# How to model the following system of ODE's in Mathematica with plots of the solution?

Consider the following

(i)

\begin{align} \frac{d S}{dt} &= -\frac{\beta S I}{N}\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} \end{align}

Where $$N=S+I$$ is the total population.

(ii)

\begin{align} \frac{dS}{dt} &= \mu N -\frac{\beta S I}{N} - \nu S\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -\nu I \end{align}

Where $$N=S+I$$ is the total population.

(iii)

\begin{align} \frac{dS}{dt} &= -\frac{\beta S I}{N} + \gamma I\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -\gamma I \end{align}

Where $$N=S+I$$ is the total population.

(iv)

\begin{align} \frac{dS}{dt} &= \mu N -\frac{\beta S I}{N} +\gamma I - \nu S\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -\gamma I -\nu I \end{align}

Where $$N=S+I$$ is the total population.

How can I model these in Mathematica with plots of the solution? What about further models with 3 or 4 equations in the system with more parameters?

A solution for the SIR model:

My attempt at a code:

ClearAll[S, P, R, b, g, tmax]
tmax = 100;
n[t_] := S[t] + P[t] + R[t]
soln = First@NDSolve[{  S'[t] == (-b/n[t])*S[t]*P[t], P'[t] == (b/n[t])*S[t]*P[t] - g P[t], R'[t] == g P[t], S[0] == 0.99, P[0] == 0.01, R[0] == 0, b == 0.07, g == 1}, {S, P, R}, {t, 0, tmax}]
Plot[{S[t] /. soln, P[t] /. soln, R[t] /. soln}, {t, 0, tmax}, PlotRange -> All, PlotLegends -> {"S", "I", "R"}, ImageSize -> Large]

• What is the issue? I have coded the first problem and I do get a solution
– user49048
Commented May 25, 2021 at 13:42
• @DiSp0sablE_H3r0 I don't know how to code it, if I get an example, I should be able code the rest(in theory anyway..)
– Math
Commented May 25, 2021 at 13:44
• Look for SIR model Commented May 25, 2021 at 13:52
• @UlrichNeumann I tried that solution but didn't get correct answer after playing with it for a bit, since my question includes $N$ with $N=S+I+R$.
– Math
Commented May 25, 2021 at 14:08
• I edited "My attempt at a code" Commented May 25, 2021 at 15:17

Here is a numerical solution of the code you presented in the comment:

tmax = 10;
b = 0.07;
g = 1;
n[t_] := S[t] + P[t] + R[t]
SPR = NDSolveValue[{S'[t] == (-b/n[t])*S[t]*P[t],
P'[t] == (b/n[t])*S[t]*P[t] - g  P[t], R'[t] == g P[t],
S[0] == 0.90, P[0] == 0.1, R[0] == 0.}, {S, P, R}, {t, 0, tmax}] ;
Plot[{Through[SPR[t]], Total@Through[SPR[t]]} // Evaluate, {t, 0,tmax}, PlotStyle -> {Automatic, Automatic, Automatic, Dashed}]


• +1 from me. I always like the solutions you suggest. This is not an exception. However, I do believe that the problem should be fully stated originally and not as a comment. This is not a comment about you, rather about the usefulness of the post itself
– user49048
Commented May 25, 2021 at 14:40
• @DiSp0sablE_H3r0 Thanks, my full agreement! Commented May 25, 2021 at 14:43
• @UlrichNeumann Can you get something similar to this? I'll attach the image in the question.
– Math
Commented May 25, 2021 at 15:02

Hopefully, after this example, you should be able to code the rest.

I am assuming that $$\beta$$ is just a constant.

Then, we begin by defining n[t_] := s[t] + i[t] as the problem suggests.

We write down the two equations:

eqn1 = D[s[t], t] + beta (s[t] i[t])/n[t];
eqn2 = D[i[t], t] - beta (s[t] i[t])/n[t];


Now, we brute-force DSolve them

DSolve[{eqn1 == 0, eqn2 == 0}, {s[t], i[t]}, t] // Flatten


A sanity check: we define the functions to be given by the answer of DSolve

s[t_] := C[1] - (E^(beta t) C[1])/(E^(beta t) - E^(C[1] C[2]))
i[t_] := (E^(beta t) C[1])/(E^(beta t) - E^(C[1] C[2]))


and then we run the two equations. If we are correct we should get two nice zeroes

eqn1 // Factor
eqn2 // Factor


which we do.

In order to plot the solutions:

we have to assign numerical values to the undetermined constants.

plot1 = Plot[{s[t]} /. C[1] -> 2 /. C[2] -> 3 /. beta -> 5, {t, 0,
5}]
plot2 = Plot[{i[t]} /. C[1] -> 2 /. C[2] -> 3 /. beta -> 5, {t, 0,
5}]
plot3 = Plot[{n[t]} /. C[1] -> 2 /. C[2] -> 3 /. beta -> 5, {t, 0, 5}]


And you get three plots.

You can combine them

Show[plot1, plot2, plot3]


• Why doesn't it let me edit the question? it says "formatting error". this happened when I asked the question which is why I had to put the equations in code format
– Math
Commented May 25, 2021 at 14:12
• I don't know what to tell you about that. Sorry. Did you try the above? Is this example clear?
– user49048
Commented May 25, 2021 at 14:13
• I tried it but didn't get a plot. also could you edit this code in the question: ClearAll[S, P, R, b, g, tmax] tmax = 100; soln = First@NDSolve[{ n[t] == S[t] + P[t] + R[t] S'[t] == (-b/n)*S[t]*P[t], P'[t] == (b/n)*S[t]*P[t] - gP[t], R'[t] == gP[t], S[0] == 0.99, P[0] == 0.01, R[0] == 0, b == 0.07, g == 1}, {S, P, R}, {t, 0, tmax}] Plot[{S[t] /. soln, P[t] /. soln, R[t] /. soln}, {t, 0, tmax}, PlotRange -> All, PlotLegends -> {"S", "I", "R"}, ImageSize -> Large]
– Math
Commented May 25, 2021 at 14:15
• Can you check again? I have shown how to plot
– user49048
Commented May 25, 2021 at 14:23
• With all due respect, you have to realize that you have not given any initial or boundary conditions to specify anything. I obtained the solutions for the problem that YOU wrote down. Suggesting that these plots do not make sense is wrong. These are the plots of the solutions for the equations of the OP. I am not a magician to guess all the missing bits of information.
– user49048
Commented May 25, 2021 at 14:37