model <- function(t, state, parms) { with(as.list(c(state,parms)), { dN1 <- r1*N1*(1-N1/K) dN2 <- r2*N2*(1-N2/K) return(list(c(dN1, dN2))) }) } p <- c(r1=1,r2=0.1,K=1) s <- c(N1=1,N2=1) run(after="state<-state+rnorm(1,mean=0,sd=0.2);if(state[1]<0)state[1]<-0") run(after="parms[\"K\"]<-rnorm(1,mean=1,sd=0.2)")