library(deSolve)
spirs <- function(t, x, parms)  {
  with(as.list(c(parms, x)), {
    dS=L-(1-u)*b*S*P+th*R-nu*S
    dI=(1-u)*b*S*P-(nu+sg+v+m+a-ph)*I
    dP= m*I-g1*P
    dR=(v+a)*I-(nu+th)*R
    output <- c(dS, dI,dP, dR)
    list(output)
  })
}
#the Initial values
#start<-c(S=0.8,P=0.08,I=0.12,R=0)
start<-c(S=1000,P=100,I=150,R=0)

## The parameters 
parms <- c(L=0.05,b=0.002,ph=0.01,a=0.01,m=0.01,nu=0.001,th=0.005,u=0, v=0,sg=0.01, g1=0.03)

## vector of timesteps
times <- seq(0,60,1) ###odel run in 52 weeks*20 years
run_d<-ode(times=times, y=start, func=spirs,parms=parms)
summary(run_d)

matplot(run_d[,2:5], type="l")
plot(run_d[,2], lwd=3, col="black", ylim=c(0,1000), type="l", main="SIR Model", ylab="Population")
lines(run_d[,3],  lwd=3,col="red")
lines(run_d[,4], lwd=3, col="blue")
lines(run_d[,5], lwd=3, col="green")
legend("topright",legend=c("S", "P","I", "R"), col=c("black","red", "blue", "green"), lty=c(1,1,1,1))

###R0
L=0.05;b=0.002;ph=0.01;a=0.01;m=0.01;nu=0.001;th=0.005;u=0; v=0;sg=0.01; g1=0.03
Ro=(1-u)*b*m*L/((g1*nu)*(m+v+a+sg+nu-ph))
Ro