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)