model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dN <- N*(b*(N/(h+N))/(1+N/k) - d)
    return(list(dN))  
  }) 
}  

p <- c(b=1,k=1,h=2,d=0.1)
s <- c(N=1)
newton()
run(state=c(N=0.35),legend=FALSE)
run(state=c(N=0.25),add=TRUE,legend=FALSE)
lines(c(0,100),c(0.3,0.3))
with(as.list(p),curve(b*(N/(h+N))/(1+N/k),from=0,to=10,xname="N",ylab="bf(N)g(N)",col="red",lwd=2))
lines(c(0,10),c(p["d"],p["d"]))

# Make Figures
size<-5#inch

pdf("alleeA.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(state=c(N=0.35),legend=FALSE)
run(state=c(N=0.25),add=TRUE,legend=FALSE)
lines(c(0,100),c(0.3,0.3))
dev.off()

pdf("alleeB.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")
with(as.list(p),curve(b*(N/(h+N))/(1+N/k),from=0,to=10,xname="N",ylab="bf(N)g(N)",col="red",lwd=2))
lines(c(0,10),c(p["d"],p["d"]))
dev.off()

pdf("alleeC.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")
plot(1,1,type="n",xlim=c(0,7),ylim=c(-0.025,0.17),xlab="N",ylab="f(N)")
with(as.list(p),curve(N*(b*(N/(h+N))/(1+N/k)-d),from=0,to=10,xname="N",ylab="bf(N)g(N)",col="red",lwd=2,add=TRUE,n=201))
lines(c(0,7),c(0,0))
dev.off()