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

p <- c(r=1,K=1,a=1,c=1,delta=0.4,n=1,e=0)
s <- c(R=1,N=0.01)
plane()
newton(c(R=0.4,N=0.6),plot=TRUE)
newton(c(R=1,N=0),plot=TRUE)
plane()
run(50,0.1,traject=TRUE)
run(50,0.1)

# Make Figures
# size<-5#inch
# 
# pdf("lotkaA.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()
# newton(c(R=0.4,N=0.6),plot=TRUE)
# newton(c(R=1,N=0),plot=TRUE)
# newton(c(R=0,N=0),plot=TRUE)
# dev.off()
# 
# pdf("lotkaB.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(legend=FALSE)
# run(50,0.1,traject=TRUE)
# dev.off()
# 
# pdf("lotkaC.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(50,0.1,legend=FALSE)
# dev.off()
# 
# pdf("lotkaD.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,a=1,c=1,delta=0.4,n=2,e=0)
# plane(legend=FALSE)
# run(50,0.1,traject=TRUE)
# dev.off()
# 
# pdf("lotkaE.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,a=1,c=1,delta=0.4,n=0.5,e=0)
# plane(legend=FALSE)
# run(50,0.1,traject=TRUE)
# dev.off()
# 
# pdf("lotkaF.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,a=1,c=1,delta=0.4,n=1,e=0)
# plane()
# p["K"]=0.6
# plane(add=TRUE,lwd=1)
# f <- newton(c(R=0.4,N=0.6),plot=TRUE)
# p["K"]=1
# run(50,0.1,state=f,traject=TRUE)
# dev.off()