model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR1 <- s - s*R1 - c11*N1*R1 - c21*N2*R1
    dR2 <- s - s*R2 - c12*N1*R2 - c22*N2*R2
    dN1 <- (b1*pmin(c11*R1/h11,c12*R2/h12) - d1)*N1
    dN2 <- (b2*pmin(c21*R1/h21,c22*R2/h22) - d2)*N2
    return(list(c(dR1, dR2, dN1, dN2)))  
  }) 
}  

qss <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    R1 <- s/(s + c11*N1 + c21*N2)
    R2 <- s/(s + c12*N1 + c22*N2)
    dN1 <- (b1*pmin(c11*R1/h11,c12*R2/h12) - d1)*N1
    dN2 <- (b2*pmin(c21*R1/h21,c22*R2/h22) - d2)*N2
    return(list(c(dN1, dN2)))  
  }) 
}  

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

p <- c(b1=0.5,b2=0.5,d1=0.1,d2=0.1,c11=0.5,c22=0.5,c12=0.25,c21=0.25,h11=1,h22=1,h12=0.1,h21=0.1,s=1)
plane(show=c("N1","N2"),zero=F);f<-run(traject=T)
newton(f,jacobian=T,rtol=0,atol=0)
plane(x="N1",y="N2",xmax=6,ymax=6,state=c(N1=0.1,N2=0.1),odes=qss,eps=-0.01)
run(state=c(N1=0.1,N2=0.1),odes=qss,traject=T)

p <- c(b1=0.5,b2=0.5,d1=0.1,d2=0.1,c11=0.5,c22=0.5,c12=0.25,c21=0.25,h11=0.5,h22=0.5,h12=0.5,h21=0.5,s=1)
plane(show=c("N1","N2"),zero=F);f<-run(400,traject=T)
newton(f,jacobian=T,rtol=0,atol=0)
plane(x="N1",y="N2",xmax=6,ymax=6,state=c(N1=0.1,N2=0.1),odes=qss,eps=-0.01,vector=T)
run(400,state=c(N1=0.1,N2=0.1),odes=qss,traject=T)