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

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


# Make Figures
# size<-5#inch
# 
# pdf("chemoMonodA.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=40)
# newton(c(R=1,N=0),plot=TRUE)
# newton(c(R=0.5,N=10),plot=TRUE)
# dev.off()
# 
# pdf("chemoMonodB.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=40,legend=FALSE)
# run(100,traject=TRUE)
# dev.off()
# 
# pdf("chemoMonodC.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,ymax=8,legend=FALSE)
# dev.off()
# 
# pdf("chemoMonodD.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")
# with(as.list(p),curve(s-w*R,from=0,to=1,xname="R",col="red"))
# for (N in seq(10)) 
# with(as.list(p),curve(a*R*N/(h+R),xname="R",col="blue",add=TRUE))
# dev.off()
# 
# pdf("chemoMonodE.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=40,show="R")
# dev.off()