How do I get outputs at fixed steps?

Question

I have a system of differential equations, there a bit complex but the closest analogy I can come to is for a catalized checmical reaction where the catalist decomposes over time, e.g.

A + B -> C + B  (rate k1)
B -> D          (rate k2)

So we have dAdt(A,B) = -k1*A*B and dBbt(a,B) = -k2*B

I can model this with rk4 and a fixed set of time points... but I want to run this until say B < 1e-8

I think this can be done by using root functions in the general ode package but that also only takes a times parameter which forces me to pick in advance the times I want the "concentrations of A and B" for.

Note: my real equations are more complex and actually quite involved to compute - they have nothing to do with chemistry other than the values are all positive numbers and once one is as close to zero as makes no odds, the state of the system stops changing.

I have a hand-rolled RK4 implementation that does what I want, but code I write is code I need to test, I know the deSolve package is well tested so I'm just looking for a way to get the outputs at fixed step sizes until a "root" function returns true

Update

Here's an example of solving my problem where I know the times I want the answer at:

kernel <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    dA <- -k1*A*B
    dB <- -k2*B
    list(c(dA,dB))
  })
}
params <- c(k1=0.03, k2=0.005)
y0 <- c(A=1.3, B=0.5)

# here I need to pick the times!
times <- seq(0,500, by=1)
out <- rk4(y0, times, kernel, params)

What I want is something like "Change the last two lines of your code to this"

out <- ___name_of_function___(
  y0, 
  initial_time=0, 
  delta_time=1, 
  kernel, 
  params, 
  rootfun=function(t,y,p){y.B<1e-8}
)

Where __name_of_function__ is hopefully some function in the deSolve package or a related helper package

Answer 1

Example 2 in https://www.rdocumentation.org/packages/deSolve/versions/1.28/topics/lsodar shows that the general approach hints that for two calls to the solver, we can have the second call evaluate all the required points, at the expense of no control over the number of evaluations in the first call:

kernel <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    dA <- -k1*A*B
    dB <- -k2*B
    list(c(dA,dB))
  })
}
params <- c(k1=0.03, k2=0.005)
y0 <- c(A=1.3, B=0.5)
rootfun <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    ifelse(B<1e-8,0,B)
  })
}

# First find where the root is
out <- ode(y0, c(0,1e8), kernel, params, root=rootfun)
# Now the second rwo of `out` will be the location of the root
# So just create the times up to that point
times <- seq(0,round(out[2,'time']), by=1)
out <- ode(y0, times, kernel, params, root = rootfun)

NOTE: The number of evaluations of the function during to root finding is limited by maxsteps which will default to 500