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

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- R*(b*R/(h+R) - d*(1+(R/k)^2))
    return(list(dR))  
  }) 
}

p <- c(b=1,k=1,h=2,d=0.15)
s <- c(R=0.1)
run(odes=model1)
run(odes=model2)

size <- 5 #inch
pdf("whalesA.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(c(s,p)),{
  curve((b/(1+R/k))*R/(h+R),xname="R",from=0,to=3.5,col="red")
  curve(d*R/R,xname="R",col="blue",add=TRUE)
})
dev.off()

pdf("whalesB.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(c(s,p)),{
  curve(b*R/(h+R),xname="R",from=0,to=1.5,col="red")
  curve(d*(1+(R/k)^2),xname="R",col="blue",add=TRUE)
})
dev.off()

pdf("whalesC.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(c(s,p)),{
  curve(R*(b/(1+R/k))*R/(h+R),xname="R",from=0,to=3.5,col="red")
  curve(R*d,xname="R",col="blue",add=TRUE)
})
dev.off()

pdf("whalesD.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(c(s,p)),{
  curve(R*b*R/(h+R),xname="R",from=0,to=1.5,col="red")
  curve(R*d*(1+(R/k)^2),xname="R",col="blue",add=TRUE)
})
dev.off()

pdf("whalesE.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(c(s,p)),curve(R*((b/(1+R/k))*R/(h+R)-d),xname="R",from=0,to=3.5,col="red"))
lines(c(0,10),c(0,0))
dev.off()

pdf("whalesF.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(c(s,p)),curve(R*(b*R/(h+R) - d*(1+(R/k)^2)),xname="R",from=0,to=1.5,col="red"))
lines(c(0,10),c(0,0))
dev.off()