How to solve a system of 4 ODEs?

I'm new to MATLAB and I'm stuck on this problem for days now. I have a model of 4 different equations and I cant figure out how to solve it. Please help me.

    betah = 1000; betam = 300; 
    gammah = 0.004; gammam = 0.06; 
    muh = 0.3; mum = 0.3; 
    alphah = 0.000353; alpham = 0.00012; 
    theta = 0.4;

These are the initial values of the constants, and the equations are given below.

    dShdt = betah - alphah*Sh*Im - gammah*Sh + theta*Ih;
    dIhdt = -muh*Ih + alphah*Sh*Im - gammah*Ih - theta*Ih;
    dSmdt = betam - alpham*Sm*Ih - gammam*Sm;
    dImdt = -mum*Im + alpham*Sm*Ih - gammam*Im; 

[Ameer Hamza]

Can you attach the variables in mathematical notation? Also, which of these variables are the constant parameters, and which are state variables?

[Shravan Nair]

The variables are Sh, Ih, Sm, Im. Rest are all constants.

Answer 1

[Star Strider]

Try this:

syms Sh(t) Ih(t) Sm(t) Im(t)  Sh0 Ih0 Sm0 Im0 T Y 
betah = 1000; betam = 300;
gammah = 0.004; gammam = 0.06;
muh = 0.3; mum = 0.3;
alphah = 0.000353; alpham = 0.00012;
theta = 0.4;


DE1 = diff(Sh) == betah - alphah*Sh*Im - gammah*Sh + theta*Ih;
DE2 = diff(Ih) == -muh*Ih + alphah*Sh*Im - gammah*Ih - theta*Ih;
DE3 = diff(Sm) == betam - alpham*Sm*Ih - gammam*Sm;
DE4 = diff(Im) == -mum*Im + alpham*Sm*Ih - gammam*Im;

[VF,Sbs] = odeToVectorField(DE1,DE2,DE3,DE4);
ODEfcn = matlabFunction(VF, 'Vars',{T,Y})

tspan = [0,50];                                                     % Choose An Appropriate Time Span Or Vector
ic = zeros(4,1)+0.1;                                                % Appropriate ?Initial Conditions
[t,y] = ode45(ODEfcn, tspan, ic);

figure
plot(t, y)
grid
legend(string(Sbs))

producing:

ODEfcn =  @(T,Y)[Y(1).*(-8.8e+1./1.25e+2)+Y(2).*Y(4).*3.53e-4;Y(1).*(2.0./5.0)-Y(2)./2.5e+2-Y(2).*Y(4).*3.53e-4+1.0e+3;Y(3).*(-3.0./5.0e+1)-Y(1).*Y(3).*1.2e-4+3.0e+2;Y(4).*(-9.0./2.5e+1)+Y(1).*Y(3).*1.2e-4];

and a plot.

[Shravan Nair]

Hey! I got a plot, thank you SO MUCH!

I was wondering if you could explain the code a bit, so that I could solve and plot ODEs on my own in the future.

What does that first line do? (One below)

syms Sh(t) Ih(t) Sm(t) Im(t)  Sh0 Ih0 Sm0 Im0 T Y

I'm assuming it defines the variables. What are the rest 6 terms for?

[Star Strider]

My pleasure!

In many situations, the easiest way to code differential equations (that can be coded) is to use the Symbolic Math Toolbox to set them up and then convert them to exercutable functions for the numeric solvers.

The syms call:

syms Sh(t) Ih(t) Sm(t) Im(t)  Sh0 Ih0 Sm0 Im0 T Y

declares the variables to be used in the subsequent code. Any terms that appear on the right-hand-side of a symbolic equation must be declared before they are used. It also declares ‘Sh(t)’ and the others as functions of ‘t’, so the diff function (in this situation) can work with them directly to produce the appropriate symbolic derivatives. It declares ‘T’ and ‘Y’ as symbolic variables that are not functions. The others such as ‘Sh0’ are not used here (I forgot to delete them), and were used in an attempt to see if dsolve could integrate the equations. Sometimes it can, however since these nonlinear equations do not have analytic solutions (some do), it cannot.

[Shravan Nair]

Aah, got it.

Then you write out the equations of the model, which is obvious. After that, what do the following lines do?

[VF,Sbs] = odeToVectorField(DE1,DE2,DE3,DE4);
ODEfcn = matlabFunction(VF, 'Vars',{T,Y})

[Star Strider]

This call to odeToVectorField:

[VF,Sbs] = odeToVectorField(DE1,DE2,DE3,DE4);

creates the system of diffferential equations as a vector of first-order differential equations, required for the numeric solvers. It also returns the ‘Sbs’ output shows the substitutions conde to create ‘VF’, so ‘Sbs(1)’ defines ‘VF(1)’ and so for the others.

This call to matlabFunction:

ODEfcn = matlabFunction(VF, 'Vars',{T,Y})

converts the vector field vector ‘VF’ to an anonymous function (here) that the numeric ODE solvers can use.

The rest of the code calls ode45 to integrate the ‘ODEfcn’ system, then plots it. The string function neatly takes the ‘Sbs’ vector and converts it to a string vector that legend can use to describe the plotted curves.

[Star Strider]

If my Answer helped you solve your problem, please Accept it!

.