model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dB <- r*B*(1 - B/K) - k*N*B
    dN <- s - d*N
    return(list(c(dB, dN)))  
  }) 
}  

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dB <- r*B*(1 - B/K) - k*N*B/(h+B)
    dN <- s - d*N
    return(list(c(dB, dN)))  
  }) 
} 

s <- c(B=0.01,N=1)
p <- c(r=1,K=1,k=0.5,d=1,s=1)
plane(ymax=2,vector=TRUE)
run(traject=TRUE)

p <- c(r=1,K=1,h=0.1,k=0.2,d=1,s=1)
plane(odes=model2,ymax=2,vector=TRUE)
run(odes=model2,traject=TRUE)
newton(c(B=0,N=1),odes=model2,plot=TRUE)
newton(c(B=0.2,N=1),odes=model2,plot=TRUE)
newton(c(B=0.8,N=1),odes=model2,plot=TRUE)

# Make figures
# size <- 5 #inch
# pdf("neutroA.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,k=0.5,d=1,s=1)
# plane(ymax=2)
# newton(c(B=0,N=1),plot=TRUE)
# newton(c(B=0.5,N=1),plot=TRUE)
# dev.off()
# 
# pdf("neutroA2.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,k=1.5,d=1,s=1)
# plane(ymax=1.5)
# newton(c(B=0,N=1),plot=TRUE)
# dev.off()
# 
# pdf("neutroB.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,k=0.2,d=1,s=1)
# plane(odes=model2,ymax=2)
# newton(c(B=0,N=1),odes=model2,plot=TRUE)
# newton(c(B=0.2,N=1),odes=model2,plot=TRUE)
# newton(c(B=0.8,N=1),odes=model2,plot=TRUE)
# dev.off()
# 
# pdf("neutroB2.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,k=0.2,d=1,s=1.75)
# plane(odes=model2,ymax=2)
# newton(c(B=0,N=1),odes=model2,plot=TRUE)
# dev.off()
# 
# pdf("neutroB3.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,k=0.2,d=1,s=0.4)
# plane(odes=model2,ymax=2)
# newton(c(B=0,N=0.5),odes=model2,plot=TRUE)
# newton(c(B=0.8,N=0.5),odes=model2,plot=TRUE)
# dev.off()