Question: Calculating and plotting Effective reproduction number in a simple example for vaccination effect in SIR model

In this exercise, we are going back to the simple SIR model, without births or deaths, to look at the effect of vaccination. The aim of this activity is to represent vaccination in a very simple way - we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination.

We are starting with the transmission and recovery parameters b = 0.4 days-1 and  c=0.1 days-1. To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.​
We are going to model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time.

I have used the following code to solve this problem but I face difficulties when I am trying to define and plot the effective reproduction number.

Any suggestions? Thank you in advance!
with(Physics, diff);


b := 0.4;
c := 0.1;
n := 10^6;
p := 0.5;
deS := diff(S(t), t) = -b*S(t)*I0(t);
deI := diff(I0(t), t) = b*S(t)*I0(t) - c*I0(t);
deR := diff(R(t), t) = c*I0(t);

F := dsolve([deS, deI, deR, S(0) = 1 - p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "Proportion of Population"], title = "SIR Model with vaccination")

#define the effective reproduction number
Reff := t -> b*1/c*S(t)*1/n;

Please Wait...