model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dS <- s - d*S - beta*S*I
    dE <- beta*S*I - (d+e)*E
    dI <- e*E - (d+delta+r)*I
    dR <- r*I - d*R
    return(list(c(dS, dE, dI, dR)))  
  }) 
}  

p <- c(s=1,d=0.01,beta=0.01,e=1,delta=0.1,r=0.1)
s <- c(S=100,E=0,I=0,R=0)
with(as.list(p), {e*beta*s/(d*(d+e)*(d+delta+r))})
run()
s["I"] <- 1
f <- run(200)
f <- newton(f)
continue(s,x="beta",y="S",log="x",xmin=1e-4,xmax=0.1,ymax=110)
continue(f,x="beta",y="S",log="x",xmin=1e-4,xmax=0.1,ymax=110,add=TRUE)

# Make Figures
size<-5#inch

pdf("seirA.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
run(100)
dev.off()

pdf("seirB.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
continue(s,x="beta",y="S",log="x",xmin=1e-4,xmax=0.1,ymax=110)
continue(f,x="beta",y="S",log="x",xmin=1e-4,xmax=0.1,ymax=110,add=TRUE)
dev.off()