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+c12*R2)/(h1+c11*R1+c12*R2) - delta1)*N1
    dN2 <- (beta2*(c21*R1+c22*R2)/(h2+c21*R1+c22*R2) - delta2)*N2
    dN3 <- (beta3*(c31*R1+c32*R2)/(h3+c31*R1+c32*R2) - delta3)*N3
    return(list(c(dR1, dR2, 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.2,beta2=0.2,beta3=0.2,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.15,c32=0.3)
pC <- c(c11=0.5,c22=0.5,c12=0.5,c21=0.5,c31=0.5,c32=0.5)  # Same diet
pH <- c(h1=0.2,h2=0.2,h3=0.2)
p <- c(pR,pN,pC,pH)

# Make Figures
size<-5#inch

pdf("parallelA.pdf",width=size,height=size)
pC <- c(c11=0.5,c22=0.5,c12=0.5,c21=0.5,c31=0.5,c32=0.5)  # Same diet
pH <- c(h1=0.2,h2=0.3,h3=0.4)
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")
colors[3]<-"red";colors[4]<-"blue";colors[5]<-"darkgreen"
plane(show=c("N1","N2","N3"),zero=F,xmax=0.8,ymax=0.8)
newton(run(1000,state=c(R1=1,R2=1,N1=0,N2=0,N3=0),timeplot=FALSE),plot=TRUE)
dev.off()

pdf("intersectA.pdf",width=size,height=size)
pC <- c(c11=0.5,c22=0.5,c12=0.24,c21=0.26,c31=0.34,c32=0.35)
pH <- c(h1=0.2,h2=0.2,h3=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")
colors[3]<-"red";colors[4]<-"blue";colors[5]<-"darkgreen"
plane(show=c("N1","N2","N3"),zero=F,xmax=0.8,ymax=0.8)
newton(run(1000,state=c(R1=1,R2=1,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
newton(run(1000,state=c(R1=1,R2=1,N1=0,N2=0.1,N3=0.1),timeplot=FALSE))
newton(run(1000,state=c(R1=1,R2=1,N1=0.1,N2=0,N3=0.1),timeplot=FALSE),plot=TRUE)
dev.off()