library(deSolve) 
sessionInfo() # check versions of R (2.11.x) and deSolve (1.7)

### the derivative function
derivs <- function(t, y, parms) {
	if( t<0) y.t.lag <- 19 # set the value of y before t = 0
	else y.t.lag <- lagvalue(t - 0.74) # calculates the value of the state variable y at t - 0.74
	
	dy <- r * y * (1-y.t.lag/m) # i.e. logistic growth with time lag
	list(dy, dy=dy) # the first component returned is used to solve for y, while the second simply keeps track of the change
	}
	
### parameters
r <- 3.5; # high r
m <- 19 # carrying capacity of 19
### initial values and times

yinit <- c(y = 19.001) # carrying capacity m is an unstable equilibrium; therefore we need to start at y not equal to 19. If we start at y=19, nothing changes.
times <- seq(-0.74, 40, by = 0.01) # we have to start time at t <= the time lag 
times <- seq(-1, 40, by = 0.01) # this works fine too.

### solve the model
yout <- dede(y = yinit, times = times, func = derivs, parms = NULL, atol = 1e-10)
## here R sees r and m. otherwise we could put r and m in a vector prms
## specify parms = prms

head(yout)

plot(yout, which = 1, type = "l", lwd = 2, main = "Lemming model", mfrow = c(1,2))
plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2)