model <- function(t, state, parms) { with(as.list(c(state,parms)), { dP <- a*P*(1 - P/K) - b*P*D/(h+P) dD <- b*P*D/(h+P) - d*D return(list(c(dP, dD))) }) } p <- c(a=0,K=0,h=0,b=0,d=0) plane(eps=-0.01) s <- c(P=0,D=0) noise <- "parms[\"K\"]<-abs(rnorm(1,1,0.1));state<-ifelse(state<0.01,0,state)" run(after=noise)