model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dN <- r*N*(1-(N/K)^m)
    return(list(dN))  
  }) 
}  

p <- c(r=1,K=1,m=1)
s <- c(N=0.01)

run(tmax=10,tstep=0.1)

# Make Figures
size<-5#inch

pdf("logistA.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.01),tmax=10,tstep=0.1,legend=FALSE,ymax=1.3)
run(state=c(N=1.25),tmax=10,tstep=0.1,legend=FALSE,add=TRUE)

p["K"] <- 1e99;colors[1] <- "blue"
run(state=c(N=0.01),tmax=6,tstep=0.1,legend=FALSE,add=TRUE,lwd=1)
p["K"] <- 1;colors[1] <- "red"
dev.off()

pdf("logistB.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(r*(1-(N/K)^m),from=0,to=1.2,xname="N",ylab="Logist",col="red",lwd=2))
lines(c(0,1.2),c(0,0))
dev.off()

pdf("logistC.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(r*(1-(N/K)^m),from=0,to=1.2,xname="N",ylab="Logist",col="red",lwd=1))
with(as.list(p),curve(r*(1-(N/K)^2),from=0,to=1.2,xname="N",ylab="Logist",col="blue",lwd=2,add=TRUE,n=201))
with(as.list(p),curve(r*(1-(N/K)^0.5),from=0,to=1.2,xname="N",ylab="Logist",col="darkgreen",lwd=2,add=TRUE,n=201))
lines(c(0,1.2),c(0,0))
dev.off()

pdf("logistD.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,1.1),ylim=c(-0.1,0.42),xlab="N",ylab="f(N)")
with(as.list(p),curve(N*r*(1-(N/K)^2),from=0,to=1.1,xname="N",col="red",lwd=2,add=TRUE,n=201))
with(as.list(p),curve(N*r*(1-(N/K)^1),from=0,to=1.1,xname="N",ylab="Logist",col="blue",lwd=1,add=TRUE,n=201))
with(as.list(p),curve(N*r*(1-(N/K)^0.5),from=0,to=1.1,xname="N",ylab="Logist",col="darkgreen",lwd=2,add=TRUE,n=201))
lines(c(0,1.2),c(0,0))
dev.off()