@mmcdara thank you so much. I am now getting what I wanted to see on my graphs.

@mmcdara below is the rcode

Rcode.txtlibrary(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

Download Rcode.txt

@mmcdara below is my maple code

Download modie-work.mw