f1 <- function(e,R) {(r*(1-R/K)*(h+R))/(a-e*r*(1-R/K))}
f2 <- function(e,R) {K*(1-a/(e*r))}

dR <- function(R,N) {r*R*(1 - R/K) - a*R*N/(h+e*N+R)}
LR <- function(R) {r*R*(1-R/K)}
FRN <- function(R,N) {a*R*N/(h+e*N+R)}
AS <- function(R) {(a/e)*R}

r <- 1; K <- 1; a <- 1; h <- 0.1; e <- 0.1

es <- c(0, 0.5, 0.75, 0.9, 0.95, 1.1)
le <- length(es) 
A <- sapply(es,f2,R)
size<-5#inch
colors <- c("red","blue","darkgreen","darkorange","darkmagenta","deeppink","gold","darkorchid","aquamarine","gray")

pdf("BeddingtonFuncA.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")
colors[1] <- "red"
curve(f1(es[1],R),xname="R",from=-0.1,to=1,ylim=c(0,2),col=colors[1],lwd=2)
for (i in seq(2,le-1)) {
  curve(f1(es[i],R),xname="R",col=colors[i],lwd=1,add=T)
  #lines(c(A[i],A[i]),c(0,3),lty=2,col=colors[i],lwd=2)
}
curve(f1(es[le],R),xname="R",col=colors[le],lwd=2,add=T)
lines(c(0,0),c(-1,5))
dev.off()

Ns <- c(0.1, 0.25, 0.5, 1, 2, 5, 10)
lN <- length(Ns)
a <- 0.5; e <- 0.8
pdf("BeddingtonFuncB.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")
curve(LR,xname="R",from=0,to=1,ylim=c(0,0.26),lwd=2)
for (i in seq(lN)) {
  curve(FRN(R,Ns[i]),xname="R",add=T)
  sol <- uniroot.all(dR,c(1e-12,1),N=Ns[i])
  if (length(sol) > 0) {
    yval <- sapply(sol,FRN,Ns[i])
    print(c(Sol=sol,Yval=yval))
    points(sol,yval,pch=16,col=colors[i],cex=2)
  }
}
curve(AS(R),xname="R",add=T,lty=2,lwd=2)
dev.off()

Ns <- c(0.1, 0.25, 0.4, 0.5, 0.6, 0.75, 1, 2, 10)
lN <- length(Ns)
a <- 0.5; e <- 0.1
pdf("BeddingtonFuncC.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")
curve(LR,xname="R",from=0,to=1,ylim=c(0,0.26),lwd=2)
for (i in seq(lN)) {
  curve(FRN(R,Ns[i]),xname="R",add=T)
  sol <- uniroot.all(dR,c(1e-12,1),N=Ns[i])
  if (length(sol) > 0) {
    yval <- sapply(sol,FRN,Ns[i])
    print(c(Sol=sol,Yval=yval))
    points(sol,yval,pch=16,col=colors[i],cex=2)
  }
}
curve(AS(R),xname="R",add=T,lty=2,lwd=2)
dev.off()