model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- r*R*(1 - R/K) - a*(R^2)*N/(h^2+R^2)
    dN <- a*(R^2)*N/(h^2+R^2) - N
    return(list(c(dR, dN)))  
  }) 
}  
dR <- function(R,N) {r*R*(1-R/K) - a*N*R^2/(h^2+R^2)}
dRL <- function(R) {r*R*(1-R/K)}
dRG <- function(R,N) {a*N*R^2/(h^2+R^2)}

r <- 1; K <- 1; a <- 1; h <- 0.1;
s <- c(R=1,N=0.01)
p <- c(r=r,K=K,h=h,a=a)

Ns <- c(0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4)
lN <- length(Ns)

size<-5#inch
colors <- c("red","blue","darkgreen","darkorange","darkmagenta","deeppink","aquamarine","gray","gold","darkorchid")

pdf("sigmoidConsa.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")
sols <- list()
curve(dRL,xname="R",from=0,to=1,ylim=c(0,0.26),lwd=2)
for (i in seq(lN)) {
  curve(dRG(R,Ns[i]),xname="R",add=T)
  sol <- uniroot.all(dR,c(1e-12,1),N=Ns[i])
  yval <- sapply(sol,dRG,Ns[i])
  print(c(Sol=sol,Yval=yval))
  if (length(sol) > 0) {
    points(sol,yval,pch=16,col=colors[i],cex=3)
    sols[[i]] <- sol
  }
}
dev.off()

pdf("sigmoidConsb.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] <- "black"
plane(xmax=1,ymax=0.4,lwd=1,show="R",legend=F)
colors[1] <- "red"
for (i in seq(length(sols))) points(sols[[i]],rep(Ns[i],length(sols[[i]])),pch=16,col=colors[i],cex=3)
dev.off()

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

r <- 1; K <- 1; a <- 1
h <- c(0.1, 0.15, 1/sqrt(28), 1/sqrt(27), 1/sqrt(26), 0.25)
lh <- length(h)


pdf("sigmoidHilla.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(R,h[1]),xname="R",ylim=c(0,0.5),col=colors[1],lwd=2)
for (i in seq(2,lh))
  curve(f1(R,h[i]),xname="R",col=colors[i],lwd=1,add=T)
lines(c(1/3,1/3),c(0,1),lty=2)
lines(c(h[1],h[1]),c(0,1),lty=2,col="red")
lines(c(0.5,0.5),c(0,1),lty=2)
dev.off()

pdf("sigmoidHillb.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(f2(R,h[1]),xname="R",from=0,to=0.5,ylim=c(-0.04,0.02),col=colors[1],lwd=2)
for (i in seq(2,lh-1))
  curve(f2(R,h[i]),xname="R",col=colors[i],lwd=1,add=T)
curve(f2(R,h[lh]),xname="R",col=colors[lh],lwd=2,add=T)
lines(c(0,0.5),c(0,0),lty=2)
lines(c(1/3,1/3),c(-0.04,0.02),lty=2)
dev.off()