model <- function(t, state, parms) { with(as.list(c(state,parms)), { b <- b0 + e*(sin(2*3.1416*(t-Delta)/365)) dA <- b*A*(1-A/k) - d*A return(list(dA)) }) } p <- c(b0=1,e=0.5,d=0.5,k=1,Delta=100) s <- c(A=1e-6) run(10*365)