I am trying to compose a MATLAB code to solve for roots using fsolve for nonlinear ODE. While trying to find the solutions for these two equations at steady state, I only obtain one solution. Is there a method of obtaining all possible steady state solutions for this problem? This is the problem and below that is the code I have so far.
ODE dx/dt=[u1-x1-72x10^9exp(-10000/x2)x1;310-x2+1.44x10^13exp(-10000/x2)x1+u2/23.9]
Given steady state uss=[0.998,0.5], find x at steady state xss=[x1,x2]
function F = root2d(x)
F(1) = 0.998 - x(1) - (((72*10^9).*exp(-10000/x(2))).*x(1));
F(2) = 310 - x(2) + (((1.44*10^13).*exp(-10000/x(2))).*x(1)) + 0.5/23.9;
fun=@root2d;
xo=[0,0];
x=fsolve(fun,xo)
x =
0.9973 310.1638
However, I saw that the solution should be xss=[0.5709,395.4047]. So I decided to use wolfram alph to solve for x, and I obtained three separate solutions for xss. Solution 1: xss=[0.00475972,508.669] Solution 2: xss=[0.5709,395.4047] Solution 3: xss=[0.9973,310.164]