model <- function(t, state, parms) { with(as.list(c(state,parms)), { dR <- r*R*(1 - R/k) - d1*R - e*R*N/(h+R) dN <- c*e*R*N/(h+R) - d2*N return(list(c(dR, dN))) }) } p <- c(r=0,k=1,h=0,e=0,c=1,d1=0,d2=0) s <- c(R=0,N=0)