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

p <- c(beta=2,r=1,s=1,d=0.01,h=10)
s <- c(S=100,I=1,R=0)
run(20)
plane(xmax=100,ymax=20)


# make figures
# size <- 5 #inch
# pdf("sifb.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")
# p <- c(beta=2,r=1,s=1,d=0.01,h=10)
# s <- c(S=100,I=1,R=0)
# plane(xmax=100,ymax=10)
# dev.off()