model <- function(t, state, parms){ with(as.list(c(state,parms)),{ N <- state[1:ndiv] dN <- rep(0,ndiv) dN[1] <- s - (p+d)*N[1] dN[2:ndiv] <- 2*p*N[1:(ndiv-1)] - (p+d)*N[2:ndiv] dS <- (1-f)*2*p*N[ndiv] - dS*S dL <- f*2*p*N[ndiv] + r*L*(1-L/K) return(list(c(dN,dS,dL))) }) } ndiv <- 25 p <- c(s=0.1,p=1/365,d=1/365,f=0.001,dS=0.1,r=1,K=100) N <- rep(0,ndiv) names(N) <- paste("N",seq(1,ndiv),sep="") s <- c(N,S=0,L=0) run(365)