The solution of SEIR model using `deSolve`

SEIR.R

library(deSolve)

SEIRmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS <- - beta*I*S
    dE <- beta*I*S - alpha*E
    dI <- alpha*E - gamma*I
    dR <- gamma*I
    return(list(c(dS, dE, dI, dR)))
  })
}

pars  <- c(alpha=0.1, beta=0.2, gamma=0.1)
yini  <- c(S=1, E=0.01, I=0, R=0)

times <- seq(0, 200, by = 1)
out   <- ode(yini, times, SEIRmod, pars)
col4 = hcl.colors(4, palette = "Set 2")
pch4 = 15:18
thin = as.integer(seq(1, nrow(out), length.out=6))
matplot(out[,1], out[,-1], type = "l", lty=1, col=col4, xlab="time", ylab = "value")
matpoints(out[thin,1], out[thin,-1], pch=pch4, col=col4)
legend("topright",legend = c("S","E", "I", "R"), lty=1, pch=pch4, col=col4)