Changing the righthand side in a callback

Hi everyone,
I am wondering how to correctly change the right hand side inside a DifferentialEquations callback.
If I don’t reinit the system runs through but gives wrong results (compared to just running new differential equations). If I reinit the integrator I get an infinite loop.
Any pointers?

mutable struct ode_switch
    original_ode
    fault_ode
    fault::Bool
end

function (pgs::ode_switch)(u, t, integrator)
    println("time $t") # To make sure the call back is really called
    if ! pgs.fault
        integrator.f = rhs(pgs.fault_ode)
        integrator.u = find_valid_ic(integrator.f, integrator.u) # This takes care of the DAE/Massmatrix aspect
        # reinit!(integrator, t0=t) # If I use reinit, I get an infinite loop with the callback being called over and over.
        # step!(integrator) # I tried to do a step in order for that not to happen but it didn't help
        pgs.fault = true
    else
        integrator.f = rhs(pgs.original_ode)
        integrator.u = find_valid_ic(integrator.f, integrator.u) # This takes care of the DAE/Massmatrix aspect
        # reinit!(integrator, t0=t) # If I use reinit, I get an infinite loop with the callback being called over and over.
        # step!(integrator) # I tried to do a step in order for that not to happen but it didn't help
        pgs.fault = false
    end
end

timespan = (0., 10.)
faulttime = (1., 2.)
pgs = ode_switch(func1,func2)
cb = FunctionCallingCallback(pgs; funcat=faultime, func_start = false)
prob = ODEProblem(func1, x0, timespan, callback=cb)
solve(prob, Rodas4(autodiff=false))

I would just use a DiscerteCallback to change a parameter that is the boolean in p, and then check for that in f. By making it integrator-local it won’t have re-use issues.

Thanks for the quick reply, I saw that as a design in the other answers. In this case it’s a bit of an awkward fit because I typically will get an ODEFunction and a fault function that changes the ODEFunction from the user. So it’s a bit of effort to build a wrapper that switches between them, thus I thought I’d try it this way.

For example, the different functions (func1, func2) could have different mass matrices, that could not be handled by a change of parameters, correct?

That probably could be, since it can be a DiffEqOperator with an update_func, but that is not well documented…

In that case I think we will stick to calling solve several times with initial conditions taken from the last simulation.

Is there a documented way to combine several solutions into one? (or the next solve on top of the existing solution object)?