model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- r*R - a*R*N - epsR*R^2
    dN <- c*a*R*N - delta*N - epsN*N^2
    return(list(c(dR, dN)))  
  }) 
}  

s <- c(R=1,N=0.01)
p <- c(r=1,a=2,c=0.5,delta=0.4,epsR=0.1,epsN=0)
plane(legend=FALSE)
newton(c(R=0.4,N=0.5),plot=TRUE)
run(50,0.1,state=c(R=0.35,N=0.35),traject=TRUE)

# Make Figures
# size<-5#inch
# 
# pdf("lotka0A.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,a=2,c=0.5,delta=0.4,epsR=0.1,epsN=0)
# plane(legend=FALSE)
# newton(c(R=0.4,N=0.5),plot=TRUE)
# run(50,0.1,state=c(R=0.35,N=0.35),traject=TRUE)
# dev.off()
# 
# pdf("lotka0B.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,a=2,c=0.5,delta=0.4,epsR=0,epsN=0)
# plane(legend=FALSE)
# run(25,0.1,state=c(R=0.45,N=0.45),traject=TRUE)
# run(25,0.1,state=c(R=0.4,N=0.4),traject=TRUE)
# run(25,0.1,state=c(R=0.35,N=0.35),traject=TRUE)
# run(25,0.1,state=c(R=0.3,N=0.3),traject=TRUE)
# run(25,0.1,state=c(R=0.25,N=0.25),traject=TRUE)
# dev.off()
# 
# pdf("lotka0C.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,a=2,c=0.5,delta=0.4,epsR=-0.1,epsN=0)
# plane(legend=FALSE)
# newton(c(R=0.4,N=0.5),plot=TRUE)
# run(40,0.1,state=c(R=0.35,N=0.35),traject=TRUE)
# dev.off()
# 
# pdf("lotka0D.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,a=2,c=0.5,delta=0.4,epsR=0,epsN=-0.1)
# plane(legend=FALSE)
# newton(c(R=0.4,N=0.5),plot=TRUE)
# run(40,0.1,state=c(R=0.35,N=0.35),traject=TRUE)
# dev.off()
# 
# pdf("lotka0E.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,a=2,c=0.5,delta=0.4,epsR=0,epsN=0.1)
# plane(legend=FALSE)
# newton(c(R=0.4,N=0.5),plot=TRUE)
# run(50,0.1,state=c(R=0.35,N=0.35),traject=TRUE)
# dev.off()
# 
# pdf("lotka0F.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,a=2,c=0.5,delta=0.4,epsR=0.2,epsN=-0.1)
# plane(legend=FALSE)
# newton(c(R=0.4,N=0.5),plot=TRUE,jacobian=TRUE)
# run(50,0.1,state=c(R=0.35,N=0.35),traject=TRUE)
# dev.off()