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

p <- c(b=1,k=1,h=0.25,d=0.25)

with(as.list(p),curve(b*N*N/(h+N),from=0,to=3,xname="N",col="red",lwd=2))
with(as.list(p),curve(d*(1+N/k)*N,xname="N",col="blue",lwd=2,add=TRUE))

# zoom in
with(as.list(p),curve(b*N*N/(h+N),from=0,to=0.25,xname="N",col="red",lwd=2))
with(as.list(p),curve(d*(1+N/k)*N,xname="N",col="blue",lwd=2,add=TRUE))

# per capita
with(as.list(p),curve(b*N/(h+N),from=0,to=4,xname="N",col="red",lwd=2))
with(as.list(p),curve(d*(1+N/k),xname="N",col="blue",lwd=2,add=TRUE))

run(state=c(N=1))
run(state=c(N=0.01),add=TRUE)


curve(b*N*N/(h+N))
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("biofilmA.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),from=0,to=4,xname="N",col="red",lwd=2))
with(as.list(p),curve(d*(1+N/k),xname="N",col="blue",lwd=2,add=TRUE))
dev.off()

pdf("biofilmB.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*N/(h+N),from=0,to=0.25,xname="N",col="red",lwd=2))
with(as.list(p),curve(d*(1+N/k)*N,xname="N",col="blue",lwd=2,add=TRUE))
dev.off()

pdf("biofilmC.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*N/(h+N),from=0,to=3,xname="N",col="red",lwd=2))
with(as.list(p),curve(d*(1+N/k)*N,xname="N",col="blue",lwd=2,add=TRUE))
dev.off()


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()