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

p <- c(r=1,K=1,a=0.5,c=1,delta=0.5)
s <- c(R=1,N=0.01)
amin <- with(as.list(p),delta/c); print(amin)

pdf("extinction.pdf",width=7,height=7)
par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0))
par(mfrow=c(2,2))

curve(with(as.list(p),delta/(c*x)),xlab="a",from=amin,to=16,ylab="R",lwd=2,col="red",main="(a)",font.main=1)
curve(with(as.list(p),(r/x)*(1-delta/(x*K))),xlab="a",from=amin,to=16,ylab="N",lwd=2,col="blue",main="(b)",font.main=1)

a <- c(0.6, 0.7, 0.8, 0.9, 1, 1.5, 2, 4, 8, 16)
n <- length(a)
org_colors <- c("gray","red","blue","darkgreen","darkorange","darkmagenta","gold","darkorchid","aquamarine","deeppink","black")
colors[1] <- org_colors[1]
colors[2] <- org_colors[1]
plane(npixels=300,xmin=-0.01,ymin=-0.01,ymax=0.6,legend=FALSE,main="(c)")
s[1] <- with(as.list(p),delta/(c*a))
s[2] <- with(as.list(p),(r/a)*(1-delta/(a*K)))
print(s)
newton(s,plot=TRUE)
for (i in seq(n)) {
  p["a"] <- a[i]
  colors[1] <- org_colors[i+1]
  colors[2] <- org_colors[i+1]
  plane(npixels=300,add=TRUE)
  s[1] <- with(as.list(p),delta/(c*a))
  s[2] <- with(as.list(p),(r/a)*(1-delta/(a*K)))
  print(s)
  newton(s,plot=TRUE)
}

p <- c(p,s=1,w=1)
curve(with(as.list(p),s*c/delta-w/x),xlab="a",from=amin,to=16,ylab="N",lwd=2,col="blue",main="(d)",font.main=1)

dev.off()