DDE events

## =============================================================================
## R-Code to solve example 7.6 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.
##
## A delay differential equation model
## Predator-prey with harvesting  - takes a while
## implemented by Karline Soetaert
## -----------------------------------------------------------------------------

library(deSolve)
LVdede <- function(t, y, p) {
 if (t > tau1) Lag1 <- lagvalue(t - tau1) else Lag1 <- yini
 if (t > tau2) Lag2 <- lagvalue(t - tau2) else Lag2 <- yini

 dy1 <- r * y[1] *(1 - Lag1[1]/K) - a*y[1]*y[2]
 dy2 <- a * b * Lag2[1]*Lag2[2] - d*y[2]

 list(c(dy1, dy2))
}

rootfun <- function(t, y, p)
  return(y[1] - Ycrit)

eventfun <- function(t, y, p)
  return (c(y[1] * 0.7, y[2]))

r <- 1; K <- 1; a <- 2; b <- 1; d <- 1; Ycrit <- 1.2*d/(a*b)
tau1 <- 0.2; tau2 <- 0.2

yini <- c(y1 = 0.2, y2 = 0.1)
times <- seq(from = 0, to = 200, by = 0.01)
yout <- dede(func = LVdede, y = yini, times = times,
           parms = 0, rootfun = rootfun,
           events = list(func = eventfun, root = TRUE))
attributes(yout)$troot [1:10]

plot(yout[,-1], type = "l")