The result of fsolve is wrong

I have system of linear equations in three unknowns.

1. x(1)-1+(1-x(3))^(M1-1)=0
2. x(2)-1+(1-x(3))^(M2-1)=0
3. x(3)-(1+x(1)^(R1+1)+I+J)*e*Pi=0

However, the results in my Matlab is wrong.

How can I do to optimize this code? Or what's wrong in my codes?

Following is my Matlab code:

function T=Analysis0828(R1,R2,W1,W2,M,a1,a2,e)
x = sym('x',[1 2 3]);

%Equations  
M1=(M/(a1/a1+a2))
M2=(M/(a2/a1+a2))

G=((W1+1)/2)*(x(1)*(1-x(1)^(R1)))/(1-x(1));
H=(x(1)^(R1+1))*((W2+1)/2)*x(2)*(1-x(2)^(R2))/(1-x(2));
Pi=(1+e+(x(1)^(R1+1))*e+(G+H)*e)^(-1);

I=x(1)*(1-x(1)^(R1))/(1-x(1));
J=x(1)*(1-x(2)^(R2))/(1-x(2));


my_eqns = [x(1)-1+(1-x(3))^(M1-1); 
           x(2)-1+(1-x(3))^(M2-1);  
           x(3)-(1+x(1)^(R1+1)+I+J)*e*Pi];  
  g = matlabFunction(my_eqns);
    F = my_vectorize(g);
x0 = [0.4,0.5,0.01];
s = fsolve(F,x0);
T=(M1*M1*(1-s(1))*s(3)+M2*M2*(1-s(2))*s(3))/M   

my_vectorize funtion code is as follows

 function f_vec = my_vectorize(fun)
    f_vec = @fun_vectorize;
    function out = fun_vectorize(x)
        x   = num2cell(x);
        out = fun(x{:});
    end
 end

When inputs is (3, 5, 5, 2, 10, 1, 1, 1)

And my results are

s =

    0.9557    0.5413    0.9557
    
    

[Walter Roberson]

my_vectorize is missing.

Also, you have T as a return value but you do not set T.

[FANG YI LIN]

Sorry, I have edited.

[Walter Roberson]

x = sym('x',[1 2 3]);

Why are you creating a 1 by 2 by 3 array of x values when you only use three of them?

[FANG YI LIN]

Thank you so so much, you point out my blind spot! (^・ω・^ )

[FANG YI LIN]

I change to the following code

x = sym('x',[1 3]);

Answer 1

[Walter Roberson]

We do not know what your my_vectorize does; in particular we do not know how it turns the function of three inputs into a function of one input. The order you arrange the inputs is important, in order for us be able to calculate which of the 37 complex-valued solutions is "closest" to [0.4,0.5,0.01]

Have you considered using the 'vars' option to matlabFunction instead of using your own custom vectorization function?

If we assume that the order for [0.4,0.5,0.01] is x1_1, x2_1, x1_2 (the first three entries of the 1 x 2 x 3 matrix of x symbols) then the one of the 37 solutions that has the smallest norm to [0.4, 0.5, 0.01] is

However, your use of the initial conditions [0.4,0.5,0.01] would imply that you want the solution that is "closest" to [0.4,0.5,0.01], and the closest is

x1_1 = 0.95573206349551101631165448445161
x2_1 = 0.95573206349551101631165448445161
x1_2 = 0.54130678130277676362675222900352

which turns out to be a re-arrangement of what you said you got as a solution. Perhaps your vectorization uses a different order.

When given the wrong order of initial conditions, fsolve() will wander off into complex values; when given the correct order of initial conditions that matches the variable in the x1_2 position being 0.01, then fsolve() will produce the numeric values 0.955732064923563 0.955732064923563 0.541306780698535

So I think that fsolve() is working -- but you have to be more careful about how you vectorize and about the order of initial conditions.