Terminate ODE solver involving root function and event function

I have taken this example from documentation for demonstration purposes. In the following example when the value of y reaches to 0.1 a random value is added. I want to terminate the solver if the y value is greater than 0.8.

One possible solution is to generate a random value in eventfun such that y is always less than 0.8. Is there any other possible solution to terminate the solver? This would be helpful in my complicated model.

## =======================================================================
## Example 3:
##   using lsodar to trigger an event
## =======================================================================

## a state variable is decaying at a first-order rate. 
## when it reaches the value 0.1, a random amount is added.
library("deSolve")

derivfun <- function (t,y,parms)
  list (-0.05 * y)

rootfun <- function (t,y,parms)
  return(y - 0.1) 

eventfun <- function(t,y,parms)
  return(y + runif(1))  


yini <- 0.5
times <- 0:400

out <- lsodar(func=derivfun, y = yini, times=times, 
              rootfunc = rootfun, events = list(func=eventfun, root = TRUE))

plot(out, type = "l", lwd = 2, main = "lsodar with event")

# }

Solution

Does the following what you want? Thanks for the clear example.

library("deSolve")

derivfun <- function (t,y,parms)
  list (-0.05 * y)

rootfun <- function (t,y,parms)
  return(c(y - 0.1, y - 0.8))

eventfun <- function(t,y,parms)
  return(y + runif(1))  


yini <- 0.5
times <- 0:400

out <- lsodar(func=derivfun, y = yini, times=times, 
              rootfunc = rootfun, 
              events = list(func=eventfun, root = TRUE, terminalroot = 2))

plot(out, type = "l", lwd = 2, main = "lsodar with event")

And here another refinement of rootfun and eventfun:

library("deSolve")

terminate <- 0.8
eps       <- 1e-6

derivfun <- function (t, y, parms)
  list (-0.05 * y)

rootfun <- function (t, y, parms)
  return(c(y - 0.1, y - terminate))

eventfun <- function(t, y, parms)
  return(min(y + runif(1), terminate + eps))  

yini <- 0.5
times <- 0:400

out <- lsodar(func=derivfun, y = yini, times=times, 
              rootfunc = rootfun, 
              events = list(func = eventfun, root = TRUE, terminalroot = 2))

plot(out, type = "l", lwd = 2, main = "lsodar with event")