model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR1 <- s1 - d1*R1 - c11*N1*R1 - c21*N2*R1
    dR2 <- s2 - d2*R2 - c12*N1*R2 - c22*N2*R2
    dN1 <- (pmin(a11*c11*R1,a12*c12*R2) - delta1)*N1
    dN2 <- (pmin(a21*c21*R1,a22*c22*R2) - delta2)*N2
    dN3 <- (pmin(a31*c31*R1,a32*c32*R2) - delta3)*N3
    return(list(c(dR1, dR2, dN1, dN2, dN3)))  
  }) 
}  

qss <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    R1 <- s1/(d1 + c11*N1 + c21*N2)
    R2 <- s2/(d2 + c12*N1 + c22*N2)
    dN1 <- (pmin(a11*c11*R1,a12*c12*R2) - delta1)*N1
    dN2 <- (pmin(a21*c21*R1,a22*c22*R2) - delta2)*N2
    dN3 <- (pmin(a31*c31*R1,a32*c32*R2) - delta3)*N3
    return(list(c(dN1, dN2, dN3)))  
  }) 
}  

s <- c(R1=1,R2=1,N1=0.1,N2=0.1,N3=0.1)

pR <- c(s1=1,s2=1,d1=1,d2=1)
pN <- c(delta1=0.25,delta2=0.26,delta3=0.25)
pC <- c(c11=1,c22=1,c12=0.5,c21=0.5,c31=0,c32=0)
pA <- c(a11=0.5,a22=0.5,a12=2,a21=2,a31=1,a32=1)
pA <- c(a11=1,a22=1,a12=1,a21=1,a31=1,a32=1)

p <- c(pR,pN,pC,pA)
initR1 <- as.numeric(p["s1"]/p["d1"]) 
initR2 <- as.numeric(p["s2"]/p["d2"])
plane(show=c("N1","N2"),xmax=1.25,ymax=1.25,zero=FALSE)
newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
run(state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),traject=TRUE)
f <- newton(run(200,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)

plane(odes=qss,state=c(N1=0,N2=0,N3=0),xmax=2,ymax=2,vector=TRUE)

# Make Figures
# size<-5#inch
# 
# pdf("tilmanMinA.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")
# pA <- c(a11=1,a22=1,a12=1,a21=1,a31=1,a32=1)
# p <- c(pR,pN,pC,pA)
# plane(show=c("N1","N2"),zero=FALSE,xmax=1.25,ymax=1.25)
# newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
# f <- newton(run(20,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
# dev.off()
# 
# pdf("tilmanMinC.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")
# plane(odes=qss,state=c(N1=0,N2=0,N3=0),xmax=2,ymax=2)
# newton(odes=qss,state=c(N1=0.5,N2=0.5,N3=0),plot=TRUE)
# dev.off()
# 
# pdf("tilmanMinB.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")
# pA <- c(a11=0.5,a22=0.5,a12=2,a21=2,a31=1,a32=1)
# p <- c(pR,pN,pC,pA)
# plane(show=c("N1","N2"),zero=FALSE,xmax=1.25,ymax=1.25)
# newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
# f <- newton(run(20,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
# dev.off()
# 
# pdf("tilmanMinD.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")
# plane(odes=qss,state=c(N1=0,N2=0,N3=0),xmax=2,ymax=2)
# newton(odes=qss,state=c(N1=0.5,N2=0.5,N3=0),plot=TRUE)
# dev.off()