model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- (r*(1 - (R/K)) - b*N/(1+R/hR+N/hR))*R
    dN <- (b*R/(1+R/hR+N/hR) - d - c*M/(1+N/hN+M/hN))*N
    dM <- (c*N/(1+N/hN+M/hN) - e)*M
    return(list(c(dR, dN, dM)))  
  }) 
}  

p <- c(r=1,K=1,b=1,c=1,d=0.4,e=0.4,hR=100,hN=100)     # Mass action
p <- c(r=1,K=1,b=10,c=10,d=0.4,e=0.4,hR=0.1,hN=0.1)   # Beddington
s <- c(R=1,N=0.01,M=0.01)
run()
plane()
s1 <- c(R=1,N=0,M=0)
s2 <- newton(run(state=c(R=1,N=0.01,M=0)))
s3 <- newton(run(state=c(R=1,N=0.01,M=0.01)))
continue(s1,x="K",y="R",positive=TRUE)
continue(s2,x="K",y="R",positive=TRUE,add=TRUE)
continue(s3,x="K",y="R",positive=TRUE,add=TRUE)

continue(s2,x="K",y="N",positive=TRUE)
continue(s3,x="K",y="N",positive=TRUE,add=TRUE)

continue(s3,x="K",y="M",positive=TRUE)

# Make Figures
size<-5#inch

p <- c(r=1,K=1,b=1,c=1,d=0.4,e=0.4,hR=100,hN=100)     # Mass action
s1 <- c(R=1,N=0,M=0)
s2 <- newton(run(state=c(R=1,N=0.01,M=0)))
s3 <- newton(run(state=c(R=1,N=0.01,M=0.01)))

pdf("chainA.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(s1,x="K",y="R",positive=TRUE)
continue(s2,x="K",y="R",positive=TRUE,add=TRUE)
continue(s3,x="K",y="R",positive=TRUE,add=TRUE)
dev.off()

pdf("chainB.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(s2,x="K",y="N",positive=TRUE)
continue(s3,x="K",y="N",positive=TRUE,add=TRUE)
dev.off()

pdf("chainC.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(s3,x="K",y="M",positive=TRUE)
dev.off()

p <- c(r=1,K=1,b=10,c=10,d=0.6,e=0.4,hR=0.1,hN=0.1)   # Beddington
s1 <- c(R=1,N=0,M=0)
s2 <- newton(run(state=c(R=1,N=0.01,M=0)))
s3 <- newton(run(state=c(R=1,N=0.01,M=0.01)))

pdf("chainD.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(s1,x="K",y="R",positive=TRUE)
continue(s2,x="K",y="R",positive=TRUE,add=TRUE)
continue(s3,x="K",y="R",positive=TRUE,add=TRUE)
dev.off()

pdf("chainE.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(s2,x="K",y="N",positive=TRUE)
continue(s3,x="K",y="N",positive=TRUE,add=TRUE)
dev.off()

pdf("chainF.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(s3,x="K",y="M",positive=TRUE)
dev.off()