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

p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.4)
s <- c(R=1,N=0.01)

p["a"] <- 0.5; p["h"] <- 0.3; plane(); run(traject=TRUE)
p["a"] <- 0.5; p["h"] <- 0.15; plane(); run(traject=TRUE)
p["a"] <- 0.75; p["h"] <- 1; plane(); run(traject=TRUE)

p["a"] <- 0.5; p["h"] <- 0.15; plane(,portrait=TRUE); run(traject=TRUE)

# Make Figures
# size<-5#inch
# 
# pdf("ratioA.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["a"] <- 0.75; p["h"] <- 1; plane(); run(traject=TRUE)
# newton(c(R=0.001,N=0.001),plot=TRUE)
# newton(c(R=0.7,N=0.6),plot=TRUE)
# dev.off()
# 
# pdf("ratioB.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["a"] <- 0.5; p["h"] <- 0.3; plane(); run(traject=TRUE)
# newton(c(R=0.001,N=0.001),plot=TRUE)
# newton(c(R=0.7,N=0.6),plot=TRUE)
# dev.off()
# 
# pdf("ratioC.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["a"] <- 0.5; p["h"] <- 0.15; plane(portrait=TRUE)
# newton(c(R=0.001,N=0.001),plot=TRUE)
# newton(c(R=0.3,N=0.6),plot=TRUE)
# dev.off()