model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- r*R*(1 - R/K) - a*R*B
    dB <- i*(Btot - B)*R - a*R*B - e*B
    return(list(c(dR, dB)))  
  }) 
}  

p <- c(r=1,K=1,a=1,i=1,Btot=2,e=0.1)
s <- c(R=1,B=0.01)
plane(xmin=-0.01,ymin=-0.01,xmax=1,ymax=1,vector=TRUE)
newton(state=c(R=0.2,B=0.8),plot=TRUE)
run(traject=TRUE)

# Make Figure
size<-5#inch

pdf("kingfisherA.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(xmin=-0.01,ymin=-0.01,xmax=1,ymax=1)
newton(state=c(R=0.2,B=0.8),plot=TRUE)
newton(state=c(R=0,B=0),plot=TRUE)
dev.off()

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- r*R*(1 - R/K) - a*R*B
    dB <- i*(Btot - B)*R/(h+R) - a*R*B - e*B
    return(list(c(dR, dB)))  
  }) 
}  

p <- c(r=1,K=1,a=1,i=1,h=1,Btot=2,e=0.1)
s <- c(R=1,B=0.01)
plane(odes=model2,xmin=-0.01,ymin=-0.01,xmax=1,ymax=1,vector=TRUE)
run(odes=model2,traject=TRUE)

pdf("kingfisherB.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(odes=model2,xmin=-0.01,ymin=-0.01,xmax=1,ymax=1)
newton(state=c(R=0.4,B=0.6),odes=model2,plot=TRUE)
newton(state=c(R=0,B=0),odes=model2,plot=TRUE)
dev.off()