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,eps=-0.01) # 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,eps=-0.01) dev.off()