model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR1 <- s1 - d1*R1 - (c11*N1 + c21*N2 + c31*N3)*R1
    dR2 <- s2 - d2*R2 - (c12*N1 + c22*N2 + c32*N3)*R2
    dN1 <- (beta1*(c11*R1/(h11+c11*R1))*c12*R2/(h12+c12*R2) - delta1)*N1
    dN2 <- (beta2*(c21*R1/(h21+c21*R1))*c22*R2/(h22+c22*R2) - delta2)*N2
    dN3 <- (beta3*(c31*R1/(h31+c31*R1))*c32*R2/(h32+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 + c31*N3)
    R2 <- s2/(d2 + c12*N1 + c22*N2 + c32*N3)
    dN1 <- (beta1*(c11*R1/(h11+c11*R1))*c12*R2/(h12+c12*R2) - delta1)*N1
    dN2 <- (beta2*(c21*R1/(h21+c21*R1))*c22*R2/(h22+c22*R2) - delta2)*N2
    dN3 <- (beta3*(c31*R1/(h31+c31*R1))*c32*R2/(h32+c32*R2) - delta3)*N3
    return(list(c(dN1, dN2, dN3)))  
  }) 
}  

s <- c(R1=1,R2=1,N1=0.11,N2=0.12,N3=0.13)

pR <- c(s1=1,s2=1,d1=1,d2=1)
pN <- c(beta1=0.4,beta2=0.4,beta3=0.4,delta1=0.1,delta2=0.1,delta3=0.1)
pC <- c(c11=0.5,c22=0.5,c12=0.25,c21=0.25,c31=0.4,c32=0.4)
pH <- c(h11=0.2,h22=0.2,h12=0.05,h21=0.05,h31=0.2,h32=0.2)
p <- c(pR,pN,pC,pH)
s <- c(R1=1,R2=1,N1=0.1,N2=0.1,N3=0.1)

plane(show=c("N1","N2","N3"),zero=F,xmax=1.2,ymax=1.2)
newton(c(R1=1,R2=1,N1=0,N2=0,N3=0),plot=T)
f<-run(10000,traject=T)
newton(f,jacobian=T,rtol=0,atol=0,plot=T)

plane(x="N1",y="N2",xmax=8,ymax=8,state=c(N1=0.1,N2=0.1,N3=0),odes=qss)
run(state=c(N1=0.1,N2=0.1,N3=0),odes=qss,traject=T)

with(as.list(p),h11/(c11*(beta1/delta1-1)))
with(as.list(p),h12/(c12*(beta1/delta1-1)))
with(as.list(p),h21/(c21*(beta2/delta2-1)))
with(as.list(p),h22/(c22*(beta2/delta2-1)))

# Make Figures
# colors[3]<-colors[1];colors[4]<-colors[2];colors[5]<-"darkgreen";size<-5#inch
# 
# pdf("essentialA.pdf",width=size,height=size)
# pC <- c(c11=0.5,c22=0.5,c12=0.25,c21=0.25,c31=0.4,c32=0.4)
# pH <- c(h11=0.2,h22=0.2,h12=0.05,h21=0.05,h31=0.2,h32=0.2)
# p <- c(pR,pN,pC,pH)
# par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
# plane(show=c("N1","N2","N3"),zero=F,xmax=1.2,ymax=1.2)
# newton(c(R1=1,R2=1,N1=0,N2=0,N3=0),plot=T)
# newton(c(R1=0.3,R2=0.3,N1=3,N2=3,N3=0),plot=T)
# dev.off()
# 
# pdf("essentialC.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(x="N1",y="N2",xmax=8,ymax=8,state=c(N1=0.1,N2=0.1,N3=0),odes=qss)
# newton(c(N1=3,N2=3,N3=0),odes=qss,plot=T)
# dev.off()
# 
# pdf("essentialB.pdf",width=size,height=size)
# pH <- c(h11=0.1,h22=0.1,h12=0.1,h21=0.1,h31=0.2,h32=0.2)
# p <- c(pR,pN,pC,pH)
# par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
# plane(show=c("N1","N2","N3"),zero=F,xmax=1.2,ymax=1.2)
# newton(c(R1=1,R2=1,N1=0,N2=0,N3=0),plot=T)
# newton(c(R1=0.3,R2=0.3,N1=3,N2=3,N3=0),plot=T)
# dev.off()
# 
# pdf("essentialD.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(x="N1",y="N2",xmax=9,ymax=9,state=c(N1=0.1,N2=0.1,N3=0),odes=qss)
# newton(c(N1=3,N2=3,N3=0),odes=qss,plot=T)
# dev.off()