model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- S + I
    dS <- -beta*S*I/N
    dI <- beta*S*I/N - delta*I
    return(list(c(dS,dI)))  
  }) 
}  

p <- c(beta=1.5,delta=0.5)
s <- c(S=6e9,I=1)
run(25,table=T)
lines(c(0,100),c(3e9,3e9))

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- S + I + E
    dS <- -beta*S*I/N
    dE <- beta*S*I/N - gamma*E
    dI <- gamma*E - delta*I
    return(list(c(dS,dE,dI)))  
  }) 
}  

p <- c(beta=1.5,gamma=2,delta=0.5)
s <- c(S=6e9,E=0,I=1)
run(40,table=T)
lines(c(0,100),c(3e9,3e9))