SEIR model not working

[RNewbie1]

Hello, I'm trying to make a basic SEIR model, but I keep getting the following errors:

1. Error in xy.coords(x, y) : 'x' and 'y' lengths differ
2. Error in SEIR.sol[, 5] : subscript out of bounds
3. Error in checkFunc(Func2, times, y, rho) :
   Model function must return a list

The output is exactly what I'm looking for, but only plots 3 lines out of the 4 that I need, below is my functional code:

[CODE] SEIR.dyn <- function(t, var, par) {

#variables and parameters

S <- var[1]

E <- var[2]

I <- var[3]

R <- var[4]

N <- S + E + I + R

beta <- par[1]

gamma <- par[2]

delta <- par[3]

#differential equations

dS <- (-betaSI)
dE <- (betaSI) - (deltaE)
dI <- (deltaE) - (gammaI)
dR <- (gammaI)

list(c(dS, dE, dI, dR))

list

} [/CODE]

And this is my running code:

[CODE] source ("SEIRmodel.R")
library(deSolve, lib.loc="C:/R")

beta <- 0.5
gamma <- 0.1
delta <- 0.5

SEIR.par<- c(beta, delta, gamma)
SEIR.init <- c(50,3,0,0)
SEIR.t <- seq(0,100, by=0.1)
SEIR.sol <- lsoda(SEIR.init, SEIR.t, SEIR.dyn, SEIR.par)

TIME <- SEIR.sol[,1]
S <- SEIR.sol[,2]
E <- SEIR.sol[,3]
I<- SEIR.sol[,4]
R<- SEIR.sol[,5]
N <- S+I

plot(TIME, S, type='l', xlab = 'Time', ylab= 'Number of individuals', col = 'purple')
lines(TIME, E, col = 'blue')
lines(TIME, I, col = 'red')
lines(TIME, R, col = 'green')
legend(20,80, legend = c("Susceptible", "Exposed", "Infected", "Recovered"), col=c("purple", "blue", "red", "green"), lty = 1) [/CODE]

Any help on how to correct these errors would be very helpful, thank you!

[pete]

Hi, It's much easier for us to help if you provide your example as a reproducible example, called a reprex .

After adding the missing *s to the differential equations (see below) it ran without errors for me, with all 4 lines plotting.

What version of dSolve are you running?

SEIR.dyn <- function(t, var, par) {
  #variables and parameters
  S <- var[1]
  E <- var[2]
  I <- var[3]
  R <- var[4]
  N <- S + E + I + R

  beta <- par[1]
  gamma <- par[2]
  delta <- par[3]
  
  #differential equations
  dS <- (-beta*S*I)
  dE <- (beta*S*I) - (delta*E)
  dI <- (delta*E) - (gamma*I)
  dR <- (gamma*I)
  
  list(c(dS, dE, dI, dR))
}