Advance ODEs One Adaptive Step at a Time in Python: SciPy RK45 step(), dense_output, and RHS switching

Advance ODEs one adaptive step at a time with SciPy RK45: use the step() method, dense_output for interpolation, and switch the RHS by reinitializing the solver.

Adaptive ODE solvers are excellent at picking stable step sizes on their own, but most high-level interfaces insist on integrating an entire time window before returning. If you need to advance the system by exactly one adaptive step at a time—because you want to react to the latest state or swap the right-hand side on each step—you can do this efficiently in Python by calling a low-level stepping interface rather than integrating through a whole interval and discarding most of the work.

What we want to achieve

Imagine a solver that lets you request a fixed number of adaptive steps and returns only those states, not the full trajectory over a user-specified grid. Conceptually it looks like this:

state = advance(rhs_fn, t_range, y0, step_count)

With step_count set to 1, the solver chooses the step size internally and returns just the next state (t1, y1) starting from (t0, y0). The next call would then start from that new point, possibly with a different rhs_fn.

Why the usual approach wastes work

High-level ODE drivers commonly integrate across a full interval and give you a dense trajectory or a vector of times. If you only ever keep the very first step, you’re asking the integrator to compute extra internal stages and accept steps you will never use. The goal here is to avoid that overhead by advancing exactly one accepted adaptive step at a time.

One-step, adaptive stepping with SciPy RK45

In SciPy, you can instantiate a solver object (for example RK45) and advance it by exactly one accepted step via its step method. The solver decides the step size adaptively and maintains its internal state between calls.

from scipy.integrate import RK45
def rhs(t, state):
    return -0.5 * state
time_marks = []
state_marks = []
t_stop = 10
time_marks.append(0.0)
state_marks.append([1.0])
engine = RK45(rhs, t0=time_marks[0], y0=state_marks[0], t_bound=t_stop)
while engine.status == 'running':
    err = engine.step()
    if err:
        print(err)
        break
    time_marks.append(engine.t)
    state_marks.append(engine.y)
import matplotlib.pyplot as plt
plt.plot(time_marks, state_marks)
plt.show()

If the plotted polyline looks jagged, that’s just the line segments connecting discrete points. When you need values between accepted steps, use interpolation from the solver’s dense output:

value_at_t = engine.dense_output()(0.75)

Updating the right-hand side each step

If you must change the ODE function at step boundaries, don’t swap it in-place on a running integrator. Solvers assume the function remains fixed; they are very sensitive to discontinuities and to any discontinuities in derivatives. If you are fine with discontinuous derivatives at the boundaries, create a new solver instance at each step so its internal setup runs again from the new initial condition with the new function.

from scipy.integrate import RK45
def drift_a(t, s):
    return -0.5 * s
def drift_b(t, s):
    return -0.2 * s
ts = [0.0]
ys = [[1.0]]
limit = 10.0
current_rhs = drift_a
integrator = RK45(current_rhs, t0=ts[-1], y0=ys[-1], t_bound=limit)
while integrator.status == 'running':
    current_rhs = drift_b if current_rhs is drift_a else drift_a
    integrator = RK45(current_rhs, t0=ts[-1], y0=ys[-1], t_bound=limit)
    err = integrator.step()
    if err:
        print(err)
        break
    ts.append(integrator.t)
    ys.append(integrator.y)
import matplotlib.pyplot as plt
plt.plot(ts, ys)
plt.show()

This pattern restarts the solver at each accepted step with a different right-hand side. Expect the output derivative to be discontinuous at the switching points. Also note that solvers perform internal initialization when they are created; when the problem changes, that setup needs to run again. In particular, if an initial step size is not provided, RK45 will re-estimate it during initialization.

Why this matters

Advancing one accepted step at a time avoids integrating through a full time window only to discard almost everything. It also gives you a clean way to implement switch-like behavior, where the system’s definition changes at specific boundaries. While this is useful for discontinuities, it does not convert a fundamentally continuous, smooth problem into something it’s not; if your system should have continuous derivatives, changing the function each step means you are solving a different, discontinuous problem.

Takeaways

When you need single-step, adaptive integration in Python, reach for a low-level solver object like RK45 and call step to advance exactly one accepted step. If you have to change the ODE right-hand side between steps, recreate the solver each time so it can reinitialize its internal state correctly. For smoother visualization or sampling, evaluate states between accepted steps using dense_output. With this pattern you keep the computation focused on the steps you actually use, and you remain in control of how and when the model changes at step boundaries.