How to use fsolve within multiple for loops?

I have a system of 4 non-linear equations with 4 unknowns x(1), x(2), x(3) and x(4). I am able to solve the equations using the following code, however, fsolve only gives me one possible solution. I determined analytically that there is one positive steady state but would like find the approximation to the SS numerically using loops for the initial guesses as well as for different parameter values and are struggling with the code for the for loops. Any help appreciated.

% The function accepts as input the variables 
%(x, m0 ,m1, gamma, C0, r1, r2, Uc, p, r4, r5, r6, a, b, r7), 
% where x = [x(1) x(2) x(3) x(4)] = [C U E M], 
% and returns vector F(x) as output.
%Parameters taken as input
m0=0.1; m1=1; gamma=1; C0=1; r1=0.01;r2=0.01; Uc=0.25; p=1.0; r4=1.0; 
r5=0.01; r6=1.0;a=0.1; b=0.1;r7=0.1;

% Calling fsolve
fun = @(x)SteadyStates(x,m0 ,m1, gamma, C0, r1, r2, Uc, p, r4, r5, r6, a, b, r7);


x0 = [0.5;0.5;0.2;0.2]; % initial values
x = fsolve(fun,x0)


function F = SteadyStates(x, m0 ,m1, gamma, C0, r1, r2, Uc, p, r4, r5, r6, a, b, r7)
C=x(1);
U=x(2);
E=x(3);
M=x(4);

F(1) = m0*(1+m1*gamma*x(1))/(1+gamma*x(1))*(C0-x(1))-r1*x(1);
F(2) = r1*x(1)+r2*x(2)*(Uc-x(2))-x(2)^2/(p+x(2))+r4*x(3)/(1+x(2))-r5*x(2)-r6*x(2)*x(3);
F(3) = a*r1*x(1)+x(2)^2/(p+x(2))-r4*x(3)/(1+x(2))-b*r6*x(2)*x(3);
F(4) = b*r6*x(2)*x(3) - r7*x(4);

end

Answer 1

[Walter Roberson]

If you have the Symbolic Toolbox, then solve F(1), F(2), F(4) for x(1), x(3), x(4) . substitute those values into F(3) and simplify, to get

(4*x2^4 + 4407*x2^3 + 43*x2^2 + 30*x2 - 40)/(4000*(x2^2 + x2 - 1))

so three of the variables can be be written in terms of x(2) which in turn gets solved as a quartic (degree 4)

You can get the complete detailed solution by solve(F,'MaxDegree',4) but I find those kinds of solutions to be essentially useless.

You can also use vpasolve() on the set of four equations. If you do that then the result is 5 solutions and a warning about possible numeric accuracy problems with one of them. But the first of the 5 solutions is a false solution that involves division by 0. You can demonstrate that the theoretical solution corresponding to the one it complains about accuracy for, is in fact an exact solution.

[Ronel Scheepers]

Thank you so much, vpasolve() worked. I now have to see if I can include a loop to vary the input of parameter values as the current values is just a randomly chosen set.