model <- function(t, state, parms) { with(as.list(c(state,parms)), { dA <- r*A*(1-A/K) - p*Z*A/(1+A) dZ <- p*Z*A/(1+A) - m*Z - F*Z^2/(h^2+Z^2); return(list(c(dA, dZ))) }) } p <- c(r=1,K=10,h=1,p=0.5,m=0.4,F=0.15) s <- c(A=10,Z=0.01) plane(xmax=20,ymax=10,eps=-0.01)