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

s <- c(R=1,N=0.01)
p <- c(r=1,K=1,h=0.1,a=0.5,c=1,e=0.1,delta=0.4)
plane();run(200,traject=TRUE);newton(c(R=0.6,N=0.6),plot=TRUE)

with(as.list(p),h/(c*a/delta-1))
with(as.list(p),K*(1-a/(e*r)))

# Make Figures
# size<-5#inch
# 
# pdf("beddingtonA.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(r=1,K=1,h=0.1,a=0.5,c=1,e=0.8,delta=0.2)
# plane(ymax=5)
# A <- with(as.list(p),K*(1-a/(e*r)))
# lines(c(A,A),c(0,5),lty=2)
# run(100,0.1,traject=TRUE)
# newton(c(R=0.5,N=1.5),plot=TRUE)
# newton(c(R=1,N=0),plot=TRUE)
# newton(c(R=0,N=0),plot=TRUE)
# dev.off()
# 
# pdf("beddingtonB.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(r=1,K=1,h=0.1,a=0.5,c=1,e=0.1,delta=0.4)
# with(as.list(p),K*(1-a/(e*r)))
# plane(legend=FALSE)
# run(200,traject=TRUE)
# newton(c(R=0.6,N=0.6),plot=TRUE)
# newton(c(R=1,N=0),plot=TRUE)
# newton(c(R=0,N=0),plot=TRUE)
# dev.off()
# 
# pdf("beddingtonC.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(r=1,K=1,h=0.1,a=0.5,c=1,e=0.1,delta=0.3)
# plane(legend=FALSE)
# run(100,0.1,traject=TRUE)
# newton(c(R=0.2,N=0.6),plot=TRUE)
# newton(c(R=1,N=0),plot=TRUE)
# newton(c(R=0,N=0),plot=TRUE)
# dev.off()
# 
# pdf("beddingtonD.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(r=1,K=1,h=0.1,a=0.5,c=1,e=0.1,delta=0.4)
# f <- newton(c(R=0.6,N=0.6))
# continue(f,x="K",y="N",xmin=0.1,xmax=4,ymin=-0.1,ymax=4,step=0.005,positive=TRUE)
# continue(c(R=1,N=0),x="K",y="N",xmin=0.1,xmax=4,ymin=-0.1,ymax=4,positive=TRUE,add=TRUE)
# p["K"] <- 2.3  # Start at Hopf bifurcation
# f <- newton(c(R=0.5,N=1.5))
# burnin <- 1e4; tmax <- 500
# for (k in seq(2.3,4,0.05)) {
#   p["K"] <- k
#   f <- run(burnin,state=f,timeplot=FALSE)
#   data <- run(tmax,tstep=0.2,state=f,table=TRUE,timeplot=FALSE)
#   nr <- nrow(data)
#   datamin <- data[(1+which(data$N[2:(nr-1)]<data$N[1:(nr-2)] & data$N[2:(nr-1)]<data$N[3:nr])),]
#   datamax <- data[(1+which(data$N[2:(nr-1)]>data$N[1:(nr-2)] & data$N[2:(nr-1)]>data$N[3:nr])),]
#   points(rep(k,nrow(datamin)),datamin$N,pch=".")
#   points(rep(k,nrow(datamax)),datamax$N,pch=".")
# }
# dev.off()
# 
# pdf("beddingtonE.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(r=1,K=1,h=0.1,a=0.5,c=1,e=0.1,delta=0.3)
# f <- newton(c(R=0.2,N=0.6))
# continue(f,x="e",y="N",xmin=0,xmax=0.5,ymin=-0.1,ymax=1.5,step=0.005,positive=TRUE)
# continue(c(R=1,N=0),x="e",y="N",xmin=0,xmax=0.5,ymin=-0.1,ymax=1.5,positive=TRUE,add=TRUE)
# p["e"] <- 0.1675  # Start at Hopf bifurcation
# f <- newton(c(R=0.3,N=0.7))
# burnin <- 1e4; tmax <- 500
# for (k in seq(0.1675,0,-0.005)) {
#   p["e"] <- k
#   f <- run(burnin,state=f,timeplot=FALSE)
#   data <- run(tmax,tstep=0.2,state=f,table=TRUE,timeplot=FALSE)
#   nr <- nrow(data)
#   datamin <- data[(1+which(data$N[2:(nr-1)]<data$N[1:(nr-2)] & data$N[2:(nr-1)]<data$N[3:nr])),]
#   datamax <- data[(1+which(data$N[2:(nr-1)]>data$N[1:(nr-2)] & data$N[2:(nr-1)]>data$N[3:nr])),]
#   points(rep(k,nrow(datamin)),datamin$N,pch=".")
#   points(rep(k,nrow(datamax)),datamax$N,pch=".")
# }
# dev.off()
# 
# pdf("beddingtonF.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(r=1,K=1,h=0.1,a=0.5,c=1,e=0.1,delta=0.3)
# f <- newton(c(R=0.2,N=0.6))
# continue(f,x="delta",y="N",xmin=0,xmax=0.5,ymin=-0.1,ymax=1.5,step=0.005,positive=TRUE)
# continue(c(R=1,N=0),x="delta",y="N",xmin=0,xmax=0.5,ymin=-0.1,ymax=1.5,positive=TRUE,add=TRUE)
# p["delta"] <- 0.35  # Start at Hopf bifurcation
# f <- newton(c(R=0.4,N=0.7))
# burnin <- 1e4; tmax <- 500
# for (k in seq(0.35,0.05,-0.01)) {
#   p["delta"] <- k
#   f <- run(burnin,state=f,timeplot=FALSE)
#   data <- run(tmax,tstep=0.2,state=f,table=TRUE,timeplot=FALSE)
#   nr <- nrow(data)
#   datamin <- data[(1+which(data$N[2:(nr-1)]<data$N[1:(nr-2)] & data$N[2:(nr-1)]<data$N[3:nr])),]
#   datamax <- data[(1+which(data$N[2:(nr-1)]>data$N[1:(nr-2)] & data$N[2:(nr-1)]>data$N[3:nr])),]
#   points(rep(k,nrow(datamin)),datamin$N,pch=".")
#   points(rep(k,nrow(datamax)),datamax$N,pch=".")
# }
# dev.off()