We will use the Euler method to solve the damped pendulum equation written as a system

for initial conditions , I'm going to assume the disease is spreading on a small island with a population of 10000 people and that it starts when a single infected person arrives on the island..

N = 10000;

S0=N-1;

I0=1;

beta=.0002;

gamma=1;

h = 1/100;

tfinal=20;

nsteps=tfinal/h;

t=linspace(0,tfinal,nsteps+1);

S=zeros(1,nsteps);

I = zeros(1,nsteps);

R = zeros(1,nsteps);

S(1)=S0; I(1)=I0;

We see that it is greater than one, so we expect an outbreak

R0=beta*N/gamma;

fprintf('The reproductive ratio is %0.2f.\n',R0)

for k=1:nsteps

[Sprime,Iprime,Rprime]=SIRode(S(k),I(k),R(k),beta,gamma);

S(k+1) = S(k) + h * Sprime;

I(k+1) = I(k) + h * Iprime;

R(k+1) = R(k) + h * Rprime;

end

S1=S;I1=I; % Saving to plot at end

plot(t,S,t,I,t,R);

xlabel('t'); legend('S','I','R')

plot(S,I,S(1),I(1),'o')

xlabel('S');ylabel('I');

fprintf('%0.1f%% of the population eventually caught the disease.\n',100*(R(end))/N)

fprintf('At the epidemic''s peak, %0.1f%% of the population was infected.\n',100*max(I)/N)

The epidemic now proceeds more slowly, so we have to run the simulation longer. We see that the total number of infected people, and the peak of the outbreak are significantly reduced.

beta=.000125;

tfinal=50;

nsteps=tfinal/h;

t=linspace(0,tfinal,nsteps+1);

S=zeros(size(t));

I = zeros(size(t));

S(1)=S0; I(1)=I0;

R0=beta*N/gamma;

fprintf('The reproductive ratio is %0.2f.\n',R0)

for k=1:nsteps

[Sprime,Iprime,Rprime]=SIRode(S(k),I(k),R(k),beta,gamma);

S(k+1) = S(k) + h * Sprime;

I(k+1) = I(k) + h * Iprime;

R(k+1) = R(k) + h * Rprime;

end

S2=S;I2=I;

plot(t,S,t,I,t,R);

xlabel('t'); legend('S','I','R')

plot(S,I,S(1),I(1),'o')

xlabel('S');ylabel('I');

fprintf('%0.1f%% of the population eventually caught the disease.\n',100*(R(end))/N)

fprintf('At the epidemic''s peak, %0.1f%% of the population was infected.\n',100*max(I)/N)

In this case, and the pandemic never takes off. As a result, only 10 people ever get sick.

beta=.00009;

tfinal=50;

nsteps=tfinal/h;

t=linspace(0,tfinal,nsteps+1);

S = zeros(1,nsteps);

I = zeros(1,nsteps);

R = zeros(1,nsteps);

S(1)=S0; I(1)=I0; R(1)=0;

R0=beta*N/gamma;

fprintf('The reproductive ratio is %0.2f.\n',R0)

for k=1:nsteps

[Sprime,Iprime,Rprime]=SIRode(S(k),I(k),R(k),beta,gamma);

S(k+1) = S(k) + h * Sprime;

I(k+1) = I(k) + h * Iprime;

R(k+1) = R(k) + h * Rprime;

end

plot(t,S,t,I,t,R);

xlabel('t'); legend('S','I','R')

plot(S,I,S(1),I(1),'o')

xlabel('S');ylabel('I');

fprintf('%0.1f%% of the population eventually caught the disease.\n',100*(R(end))/N)

fprintf('At the epidemic''s peak, %0.1f%% of the population was infected.\n',100*max(I)/N)

fprintf('In absolute terms, at the epidemic''s peak there were %i people infected and a total of %i people ever got infected.\n',ceil(max(I)),ceil(max(R)))

As computed above, in the first two experiments but in the third experiment, which correlates with our observation that outbreaks occur in the first two experiments but not in the third.

We plot the susceptible and infected populations in all three experiments. The result of experiment 3 barely visibile on this axis since the number of infected individuals goes to zero so fast and the the number of susceptible barely changes.

plot(S1,I1,S2,I2,S,I,'o')

xlabel('S');ylabel('I')

legend('Experiment 1','Experiment 2','Experiment 3 ')

function[Sprime,Iprime,Rprime]=SIRode(S,I,R,beta,gamma)

Sprime = -beta*S*I;

Iprime = beta*S*I - gamma*I;

Rprime = gamma*I;

end