model <- function(t, state, parms) { with(as.list(c(state,parms)), { vR <- pmax(0,R-f) #Or: vR <- ifelse(R-f < 0, 0, R-f) dR <- r*R*(1 - R/K) - a*vR*N/(h+vR) dN <- c*a*vR*N/(h+vR) - delta*N return(list(c(dR, dN))) }) } p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.4,f=0.05) s <- c(R=1,N=0.01) plane(eps=-0.01)