model <- function(t, state, parms) { with(as.list(c(state,parms)), { dW <- R*(w0 + V/(k2 + V)) - rW*W - g*V*W/(k1 + W) dV <- i + c*g*V*W/(k1 + W) - d*V - h*V return(list(c(dW, dV))) }) } p <- c(R=2,c=10,d=0.25,g=0.05,i=0.01,k1=5,k2=5,rW=0.2,w0=0.2,h=0) s <- c(W=5,V=100) plane(0,10,0.01,1000,log="y")