MATLAB: How to use fsolve within multiple for loops

for loopfsolveinitial guesses

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

Best Answer

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.
Related Question