library(deSolve)

parameters<-c(b=1.0, d=0.2, beta=0.1, gamma=0.7, D=0.3)
state<-c(X=15, Y=4, Z=0)

SIR<- function(t,state,parameters){
   with(as.list(c(state,parameters)),{
# Equations
      dX = b*(d*X+D*Y+d*Z) - beta*X*Y - d*X
      dY = beta*X*Y - (D+gamma)*Y
      dZ = gamma*Y -d*Z
      return(list(c(dX,dY,dZ)))
   })
}

times<- seq(1,100,by=1)
out<-as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(out)

par(mfrow=c(2,2), oma=c(0,0,3,0))
plot(times,out$X, type="l", main="Susceptible", xlab="time", ylab="-")
plot(times,out$Y, type="l", main="Infecté", xlab="time", ylab="-")
plot(times,out$Z, type="l", main="Guéri", xlab="time", ylab="-")