model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- s - w*R - a*R*N
    dN <- c*a*R*N - delta*N
    return(list(c(dR, dN)))  
  }) 
}  

p <- c(s=1,w=1,a=0.1,c=20,delta=1)
s <- c(R=1,N=0.01)
plane(ymax=50)
newton(c(R=1,N=0),plot=TRUE)
newton(c(R=0.5,N=10),plot=TRUE)
plane(ymax=50)
run(20,traject=TRUE)
run(20,0.1)

p["s"] <- 1
plane(ymax=30,legend=FALSE)
p["s"] <- 0.7
plane(ymax=30,add=TRUE,lwd=1)
f <- newton(c(R=0.5,N=5),plot=TRUE)
p["s"]=1
run(50,0.1,state=f,traject=TRUE)

# Make Figures
# size<-5#inch
# pdf("chemoA.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")
# plane(ymax=50)
# newton(c(R=1,N=0),plot=TRUE)
# newton(c(R=0.5,N=10),plot=TRUE)
# dev.off()
# 
# pdf("chemoB.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")
# plane(ymax=50,legend=FALSE)
# run(20,traject=TRUE)
# dev.off()
# 
# pdf("chemoC.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(20,0.1,ymax=12,legend=FALSE)
# dev.off()
# 
# pdf("chemoD.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(s=1,w=1,a=0.1,c=20,delta=1)
# plane(ymax=30)
# p["s"] <- 0.7
# plane(ymax=30,add=TRUE,lwd=1)
# f <- newton(c(R=0.5,N=5),plot=TRUE)
# p["s"]=1
# run(50,0.1,state=f,traject=TRUE)
# dev.off()