DDE

## =============================================================================
## R-Code to solve example 7.3 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.  
## Solving differential equations in R. UseR, Springer, 248 pp. 
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
## implemented by Karline Soetaert
## =============================================================================

## A delay differential equation model
library(deSolve)

mackey <- function(t, y, parms, tau) {
  tlag <- t - tau
  if (tlag <= 0)
    ylag <- 0.5
  else 
    ylag <- lagvalue(tlag)
  dy <- 0.2 * ylag * 1/(1+ylag^10) - 0.1 * y
  list(dy = dy, ylag = ylag)
}

yinit <- 0.5
times <- seq(from = 0, to = 300, by = 0.1)

yout1 <- dede(y = yinit, times = times, func = mackey, 
             parms = NULL, tau = 10)
yout2 <- dede(y = yinit, times = times, func = mackey, 
             parms = NULL, tau = 20)

plot(yout1, lwd = 2, main = "tau=10",
    ylab = "y", mfrow = c(2, 2), which = 1)
plot(yout1[,-1], type = "l", lwd = 2, xlab = "y") 
plot(yout2, lwd = 2, main = "tau=20",
    ylab = "y", mfrow = NULL, which = 1)
plot(yout2[,-1], type = "l", lwd = 2, xlab = "y")