In[]:=
Needs["RobertNachbar`EpidemiologyModeling`"]
In[]:=
eqns={
′
s
[t]Λ-μs[t]-βs[t]λ[t],​​
′
i
[t]-γi[t]-μi[t]+βs[t]λ[t],​​
′
r
[t]γi[t]-μr[t]};
In[]:=
vars={s,i,r};
In[]:=
foi={λ[t]i[t]/n};
In[]:=
lambdaSoln=First@Solve[Thread[Total@eqns/.
′
_
[t]->0,Equal],Λ]//Simplify//ReplaceAll[Total@Through[vars[t]]->n]
Out[]=
{Λnμ}
In[]:=
ngm=NextGenerationMatrix[​​eqns/.lambdaSoln,​​vars,​​t,​​foi,​​{i}​​]
Out[]=

β
γ+μ

In[]:=
brn=
ℜ
0
(Max@Abs@Eigenvalues[ngm]//Simplify[#,{β,γ,μ}∈PositiveReals]&)
Out[]=
ℜ
0

β
γ+μ
In[]:=
betaSoln=First@Solve[brn,β]
Out[]=
{β(γ+μ)
ℜ
0
}
In[]:=
params=<|
ℜ
0
->1.90,β->β,γ->0.035,μ->0.07,Λ->Λ,n->10334095,I0->1|>/.lambdaSoln/.betaSoln
Out[]=

ℜ
0
1.9,β(γ+μ)
ℜ
0
,γ0.035,μ0.07,Λnμ,n10334095,I01
In[]:=
paramsN=Normal@(params//.params)
Out[]=
{
ℜ
0
1.9,β0.1995,γ0.035,μ0.07,Λ723387.,n10334095,I01}
In[]:=
inits={s[0]==n-1,​​i[0]==1,​​r[0]==0};
In[]:=
tmax=350;
In[]:=
soln=NDSolve[{eqns,inits}/.lambdaSoln/.foi/.paramsN,vars,{t,0,tmax}];
In[]:=
fig1=Plot[Evaluate[Through[vars[t]]/.soln],{t,0,tmax},PlotLegends->{"Susceptible","Infected","Recovered"},AxesLabel->{"Time","Population"},AxesLabel->{"Time","Population"},PlotStyle->{Blue,Red,Green},PlotRange->All]
Out[]=
Susceptible
Infected
Recovered