seirds.f <- function(t, y, k) {
  S <- y[1]
  E <- y[2]
  I_R <- y[3]
  I_D <- y[4]
  R <- y[5]
  D <- y[6]
  N <- y[7]

  beta = k[1]
  delta = k[2]
  gamma_r = k[3]
  gamma_d = k[4]
  mu = k[5]
  epsilon = k[6]
  waning = k[7]

  dS <- S - beta*(I_R+I_D)/N + waning*R
  dE <- E + beta*S*(I_R+I_D)/N - delta*E + epsilon
  dI_R <- I_R + delta*(1-mu)*E - gamma_r*I_R + epsilon
  dI_D <- I_D + delta*mu*E - gamma_d*I_D + epsilon
  dR <- R + gamma_r*I_R - waning*R
  dD <- D + gamma_d*I_D
  return(as.matrix(c(dS, dE, dI_R, dI_D, dR, dD, N)))
}

# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0
seirds.params <- c(0.7316455696202532,    # beta
                   0.020253164556962026,  # delta
                   0.09113924050632911,   # gamma_r
                   0.002531645569620253,  # gamma_d
                   0.0034602076124567475, # mu
                   0.0)                   #epsilon

seirds.params <- c(1,    # beta
                   1,  # delta
                   1,   # gamma_r
                   1,  # gamma_d
                   1, # mu
                   1)                   #epsilon


S <- 62
E <- 8
I_R <- 288
I_D <- 1
R <- 36
D <- 1

N <- S+E+I_R+I_D+R+D
S <- S/N
E <- E/N
I_R <- I_R/N
I_D <- I_D/N
R <- R/N
D <- D/N

tmin <- 0
tmax <- 20

install.packages("pracma")
library(pracma)
seirds.ode.sol <- ode45(f = function(t,y){seirds.f(t,y,k=seirds.params)},
                     y = c(S, E, I_R, I_D, R, D, N),
                     t0 = tmin, tfinal = tmax)
# plot
plot.seirds <- function(sol, method){
  # - pred/prey columns melted to 1 column called "value"
  sol.melted <- melt(data.frame(time=sol$t,
                                S=sol$y[,1],
                                E=sol$y[,2],
                                I_R=sol$y[,3],
                                I_D=sol$y[,4],
                                R=sol$y[,5],
                                D=sol$y[,6],),
                     id = "time") # melt based on time
  # New column created with variable names, called “variable”
  colnames(sol.melted)[2] <- "Group" # used in legend
  g <- ggplot(data = sol.melted,
              aes(x = time, y = value, color = Group))
  g <- g + geom_point(size=2)
  g <- g + xlab("time") + ylab("Population")
  g <- g + ggtitle(paste("SEIRDS Model Using", method))
  show(g)
}

plot.seirds(seirds.ode.sol, "ODE45")