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)