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)))  
  }) 
}  

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+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(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.5,c21=0.5,c31=0.5,c32=0.5)  # Same diet
pH <- c(h1=0.2,h2=0.2,h3=0.2)
pC["c11"] <- rnorm(1,0.5,0.1)
pC["c22"] <- rnorm(1,0.5,0.1)
pC["c12"] <- rnorm(1,0.25,0.1)
pC["c21"] <- rnorm(1,0.25,0.1)
pC["c31"] <- rnorm(1,0.75/2,0.1)
pC["c32"] <- rnorm(1,0.75/2,0.1)
p <- c(pR,pN,pC,pH)

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)
newton(run(1000,state=c(R1=1,R2=1,N1=0.1,N2=0,N3=0),timeplot=FALSE),plot=TRUE)
newton(run(1000,state=c(R1=1,R2=1,N1=0,N2=0.1,N3=0),timeplot=FALSE),plot=TRUE)
newton(run(1000,state=c(R1=1,R2=1,N1=0,N2=0,N3=0.1),timeplot=FALSE),plot=TRUE)
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),plot=TRUE)
newton(run(1000,state=c(R1=1,R2=1,N1=0.1,N2=0,N3=0.1),timeplot=FALSE),plot=TRUE)

colors[1]<-"red";colors[2]<-"blue";colors[3]<-"darkgreen"
cube(x=3, y=1, z=2, theta=130,xmax=7,ymax=9,zmax=9,odes=qss,state=c(N1=0.1,N2=0.1,N3=0.1))