# Source grind.R and cube.R first!

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
     dR  <- R*(r*(1-R) - a1*N1/(h+R) - a2*N2)
     dN1 <- N1*(a1*R/(h+R) - d)
     dN2 <- N2*(a2*R - d)
    return(list(c(dR, dN1, dN2)))  
  }) 
}  

p <- c(r=1,a1=1,a2=1,h=1,d=1)
s <- c(R=0.1,N1=0.1,N2=0.1)

with(as.list(p),{curve(a1*R/(h+R),from=0,to=1,xname="R",ylim=c(0,1),col="red",lwd=2,ylab="f(R)")})
with(as.list(p),{curve(a2*R,xname="R",col="blue",lwd=2,add=TRUE)})

plane(ymax=2,zero=FALSE,show=names(s))

cube(ymax=2,zmax=1.5)
cube(shade="z",show="R",add=TRUE)
f <- run3d(tmax=1000,tstep=0.1,add=TRUE)
plotdev(theta=20,phi=25) 
run(100,0.1,state=f)