Youor description is essentially correct. However, odeToVectorField does not ‘simplify’ the system, instead it converts it to a vector of first-order differential equations that — after converting them to an anonymous function — the MATLAB numeric ODE integration functions can use.
‘My question is what is the order of the results.’
That is provided in the ‘Subs’ output of odeToVectorField.
Taking the posted code only as far as the ode45 call, and plotting the results demonstrates this —
syms S(t) P(t) c(t) z(t) A b g h v k l u d r w t Y
eqn1 = diff(S,t) == (A*exp(h*t)*b*S)^(1-g) * z^g - v*z - c - A*exp(h*t)*b*S*(k*exp(h*t*(-1))+1);
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S;
eqn3 = diff(c,t) == ( ((A*exp(h*t)*b) * ((1-g)*b*z^g*(A*exp(h*t)*S*b)^(-g)*(A*exp(h*t)+S) - 1 - k*exp(l*t*(-1)))) + (((1-b)/u) * (v - g*z^(g-1)*(A*exp(h*t)*S*b)^(1-g))) - r ) * c ;
eqn4 = diff(z,t) == ( ( (g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1) - v)/(g*(g-1)*(A*exp(h*t)*S*b)^(1-g)*z^(g-2)) ) * ( d + ((c*u*(1-w))/(P*w*(v-g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1)))) + ((A*exp(h*t)*b) * ((1-g)*z^g*(A*exp(h*t)*b*S)^(-g)*(A*exp(h*t)*b + b*S) - 1 - (k*exp(l*t*(-1))))) + (((1-b)*(v-g*(A*exp(h*t)*b*S)^(1-g)*z^(g-1)))/(u)) ) ) + ( ( (b*(A*exp(h*t)+S)*z)/(A*exp(h*t)*S) ) * ( S*A*exp(h*t)*h + (A*exp(h*t)) * ((A*exp(h*t)*b*S)^(1-g)*z^g - v*z - c - (1+(k*exp(l*t*(-1))))*A*exp(h*t)*b*S) ) );
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,g,h,k,l,r,u,v,w});
A=3;
b=0.4;
g=0.3;
h=0.3;
v=2;
k=4;
l=0.7;
u=10;
d=0.2;
r=0.5;
w=0.7;
tspan = [0 10];
y0 = [100 60 30 50];
[t,Y] = ode15s(@(t,Y) odefcn(t,A,Y,b,d,g,h,k,l,r,u,v,w) , tspan, y0);
figure
hp = plot(t,real(Y));
hold on
plot(t, imag(Y),'--')
hold off
grid
legend(hp, string(Subs))
Also, note that ode15s is more appropriate for this, since the system appears to be ‘stiff’.
.
