# Solving and plotting an SIR epidemiology model

The model is given by:

Sus'[t] = (-a) Sus[t] Inf[t] ,

Inf'[t] = a Sus[t] Inf[t] - b Inf[t] ,

Recov'[t] = b Inf[t],

Sus = 'some number' , Inf = 'some number', Recov = 'some number' ,

a = 'some number', b = 'some number'


I did not find any help on-line and the explanation in the lecture is far from clear for me.

I tried to adapt the code from this question, but failed and was asked to use the former one.

• This demonstration should answer this question: "SIR Epidemic Dynamics". – Anton Antonov Oct 23 at 14:46
• 1. "I am not looking for theoretical aspects, but for a chunk of code to plot it. " -- Please, examine that demonstration again -- the code is in the demonstration. 2. "This is the problem of this tool - not enough resources to educate myself." -- What tool are you talking about? – Anton Antonov Oct 23 at 14:57
• "I am very new to Mathematica. And I find it much less abundant with code examples than R, which I usually use." -- Mathematica is much better documented than R -- since R is often somewhat inadequate people have written and posted on the web a lot about how to deal with R's inadequacies. (If you curious about how and why that happens find and read the article "The Lisp Curse".) – Anton Antonov Oct 23 at 15:24
• "I do not quite understand where is it in the 'demonstration'?" -- Two ways to get the code. 1) Download the demonstration notebook and in that notebook press the button "Download source code". 2) In the demonstration's web page page press the button "Source". – Anton Antonov Oct 23 at 15:26
• Good luck with your studies! (And my code is pretty straightforward application following the documentation.) – Anton Antonov Oct 24 at 8:26

I more or less copied and pasted into NDSolve the formulation given in the original post and replaced "=" with "==".

ClearAll[Sus, Inf, Recov, a, b, tmax]

tmax = 20;

soln = First@NDSolve[{
Sus'[t] == (-a)*Sus[t]*Inf[t],
Inf'[t] == a*Sus[t]*Inf[t] - b*Inf[t],
Recov'[t] == b*Inf[t],
Sus == 762,
Inf == 1,
Recov == 0, a == 0.00218, b == 0.44036}, {Sus, Inf, Recov}, {t,0, tmax}]

Plot[{Sus[t] /. soln, Inf[t] /. soln, Recov[t] /. soln}, {t, 0, tmax},
PlotRange -> All, PlotLegends -> {"Sus", "Inf", "Recov"},
ImageSize -> Large] 