# How to numerically solve this system?

Our system:

\begin{align} \dot S &= \mu -\beta_1 S I_2 -\beta_2 S J-\beta_3 S A - \nu S\\[2ex] \dot {I_1} &= p \beta_1 S I_1 +q \beta_2 S J +r \beta_3 S A + \xi_1 J- b_1 I_1 \\[2ex] \dot{I_2} &=(1 - p) \beta_1 S I_2 + (1 - q) \beta_2 S J + (1 - r) \beta_3 S A + \epsilon I_1+ \xi_2 J - b_2 I_2 \\[2ex] \dot{J} &=p_1 I_2 -b_3 J \\[2ex] \dot{A} &= p_2 J -b_4 A \end{align}

With $$N= S+I_1+I_2 +J+A$$. Implying $$\dot N = \mu - \nu N -\alpha A \leq \mu-\nu N$$

I tried:

tmax = 1000;
\[Beta]1 = 0.2;
\[Beta]2 = 0.1;
\[Beta]3 = 0.1;
p = 0.8;
q = 0.9;
r = 0.7;
\[Xi]1 = 0.1;
\[Xi]2 = 0.3;
\[Epsilon] = 0.2;
b1 = 0.4;
b2 = 0.3;
b3 = 0.5;
b4 = 0.2;
p1 = 0.3;
p2 = 0.2;
\[Mu] = 0.3;
\[Nu] = 0.1;
SIIJA = NDSolveValue[{
S'[t] == \[Mu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*
J[t] - \[Beta]3*S[t]*A[t] - \[Nu]*S[t],
I1'[t] ==
p*\[Beta]1*S[t]*I2[t] + q*\[Beta]2*S[t]*J[t] +
r*\[Beta]3*S[t]*A[t] + \[Xi]1 *J[t] - b1*I1[t],
I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
J[t] + (1 - r)*\[Beta]3*S[t]*A[t] + \[Epsilon]*I1[t] + \[Xi]2*
J[t] - b2*I2[t],
J'[t] == p1*I2[t] - b3*J[t],
A'[t] == p2*J[t] - b4*A[t] ,
S[0] == 50,
I1[0] == 50,
I2[0] == 50,
A[0] == 10,
J[0] == 10},
{S, I1, I2, J, A},
{t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax},
PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True,
FrameLabel -> st /@ {"Time", "Density"},
PlotLegends ->
Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)",
"I1(t)", "I2(t)", "J(t)", "A(t)" },
LegendFunction -> Framed], {0.85, 0.65}], ImageSize -> 500]


but this doesn't seem correct..

EDIT:

Our "full" system:

    tmax = 1000;
\[Beta]1 = 0.00001;
\[Beta]2 = 0.0006;
\[Beta]3 = 0.00001;
p = 0.9;
q = 0.8;
r = 0.3;
\[Xi]1 = 0.8;
\[Xi]2 = 0.9;
\[Epsilon] = 0.0002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.55;
\[Nu] = 0.01;
SIIJA = NDSolveValue[{
S'[t] == \[Mu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*
J[t] - \[Beta]3*S[t]*A[t] - \[Nu]*S[t],
I1'[t] ==
p*\[Beta]1*S[t]*I2[t] + q*\[Beta]2*S[t]*J[t] +
r*\[Beta]3*S[t]*A[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Nu])*I1[t],
I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
J[t] + (1 - r)*\[Beta]3*S[t]*A[t] + \[Epsilon]*I1[t] + \[Xi]2*
J[t] - (p1 + \[Nu])*I2[t],
J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Nu])*J[t],
A'[t] == p2*J[t] - (\[Alpha] + \[Nu])*A[t] ,
S[0] == 50,
I1[0] == 50,
I2[0] == 50,
J[0] == 10,
A[0] == 10},
{S, I1, I2, J, A},
{t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax},
PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True,
FrameLabel -> st /@ {"Time", "Density"},
PlotLegends ->
Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)",
"I1(t)", "I2(t)", "J(t)", "A(t)" },
LegendFunction -> Framed], {0.85, 0.75}], ImageSize -> 750]


EDIT 2:

Their plots(in the paper attached) are different to the plots I got when I modelled the same system, is there a reason why?

tmax = 2000;
\[Beta]1 = 0.0001;
\[Beta]2 = 0.006;
p = 0.9;
q = 0.8;
\[Xi]1 = 0.8;
\[Xi]2 = 0.9;
\[Epsilon] = 0.002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.01;
\[Nu] = 0.55;
SIIJA = NDSolveValue[{
S'[t] == \[Nu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*J[t] - \[Mu]*
S[t],
I1'[t] ==
p*\[Beta]1*S[t]*I2[t] +
q*\[Beta]2*S[t]*J[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Mu])*I1[t],
I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
J[t] + \[Epsilon]*I1[t] + \[Xi]2*J[t] - (p1 + \[Mu])*I2[t],
J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Mu])*J[t],
A'[t] == p2*J[t] - (\[Alpha] + \[Mu])*A[t] ,
S[0] == 50,
I1[0] == 50,
I2[0] == 50,
J[0] == 20,
A[0] == 20},
{S, I1, I2, J, A},
{t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax},
PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True,
FrameLabel -> st /@ {"Time", "Density"},
PlotLegends ->
Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)",
"I1(t)", "I2(t)", "J(t)", "A(t)" },
LegendFunction -> Framed], {0.85, 0.45}], ImageSize -> 500]


Here is the reproduction number:

r0 = ((\[Beta]1 b3 (\[Epsilon] p +
b1 (1 - p)) + \[Beta]2 p1 (\[Epsilon] q +
b1 (1 - q))) (\[Nu]/\[Mu]))/(b1 b2 b3 - (\[Epsilon] \[Xi]1 +
b1 \[Xi]2) p1);


EDIT 3

tmax = 2000;
\[Beta]1 = 0.0001;
\[Beta]2 = 0.006;
p = 0.3;
q = 0.4;
\[Xi]1 = 0.001;
\[Xi]2 = 0.003;
\[Epsilon] = 0.0002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.01;
\[Nu] = 0.55;
SIIJA = NDSolveValue[{
S'[t] == \[Nu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*J[t] - \[Mu]*
S[t],
I1'[t] ==
p*\[Beta]1*S[t]*I2[t] +
q*\[Beta]2*S[t]*J[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Mu])*I1[t],
I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
J[t] + \[Epsilon]*I1[t] + \[Xi]2*J[t] - (p1 + \[Mu])*I2[t],
J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Mu])*J[t],
A'[t] == p2*J[t] - (\[Alpha] + \[Mu])*A[t] ,
S[0] == 50,
I1[0] == 50,
I2[0] == 50,
J[0] == 20,
A[0] == 20},
{S, I1, I2, J, A},
{t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax},
PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True,
FrameLabel -> st /@ {"Time", "Density"},
PlotLegends ->
Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)",
"I1(t)", "I2(t)", "J(t)", "A(t)" },
LegendFunction -> Framed], {0.85, 0.75}], ImageSize -> 500]


We know a priori that $$N(2000)$$ should be approximately $$55$$. however from the model above if we do f1[2000] + f2[2000] + f3[2000] + f4[2000] + f5[2000] it gives 51.3631. This is what I was trying to explain to Ulrich.

• Correct the typo, double I2, in SIIJA : {S, I2, I2, J, A}  Nov 29 '21 at 14:55
• @Akku14 Thanks, however $A$ isn't still showing up? Also, I don't know how to include the condition $N=S+I_1+I_2+J+A$
– Math
Nov 29 '21 at 14:57
• There is a typo in your Plot-command: Change f[5] to f5[t] Nov 29 '21 at 14:59
• @Math S'[0] == -604.5 so it's not surprising that S[t] quickly approaches zero. Nov 29 '21 at 15:46
• @UlrichNeumann How do I verify that $S+I_1+I_2+J+A$ equals to $N$ at each point in time?
– Math
Nov 29 '21 at 17:02