I am using OrdinaryDiffEq to solve a set of ordinary differential equations. During the time integration I would like to extend the solution vector and change the problem. For that I use a callback of type DiscreteCallback and reinitialize the integrator with the new problem:

using OrdinaryDiffEq

function modify_integrator!(integrator)
    # Define the extended ODE.
    function growth!(du, u, p, t)
        println("growth!")
        du[1] = 0.1*u[1]
        du[2] = 0.2*u[2]
        du[3] = 0.3*u[3]
    end
    u0 = [integrator.u[1]; integrator.u[2]; 1.0]
    tspan = (integrator.t, 10.0)

    # Define the modified problem.
    prob = ODEProblem(growth!, u0, tspan, callback=callbacks)

    # Set up the integrator for the modified problem.
    iter_old = integrator.iter
    naccept_old = integrator.stats.naccept
    integrator = init(prob, integrator.alg, dt=integrator.dt, callback=integrator.opts.callback)
    integrator.iter = iter_old
    integrator.stats.naccept = naccept_old
end


# Define the ODE.
function growth!(du, u, p, t)
    du[1] = 0.1*u[1]
    du[2] = 0.2*u[2]
end

# Define the initial conditions and time span.
u0 = [1.0; 1.0]
tspan = (0.0, 10.0)

# Define the callbacks.
condition(u, t, integrator) = t > 9
cb = DiscreteCallback(condition, modify_integrator!, save_positions=(false,false))
callbacks = CallbackSet(cb)

prob = ODEProblem(growth!, u0, tspan, callback=callbacks)
sol = solve(prob, Tsit5())

The integrator starts with a u-vector with 2 components and the end result should be a vector with 3 components. But my end result is as if the modified integrator got ignored entirely and the original problem is being solved. In my debugging attempts I could verify that the modified growth! function is indeed called, but with no effect on the end solution. It is as if ODEProblem retained the original problem somehow.

The function reinit! seems like it could almost help me out here, but I don’t see how to give it a modified problem.