statefully evolving an auxiliary variable

Given any vector field f(t, y, args), I want to keep track of an auxiliary variable w that depends and evolves with y but non-differentiably in the sense that there is no vector field for w, it's just a state that gets updated every step call. For example, at every step of the integrator I would ideally like to be able to do something like the following pseudocode:

def update_w(t: float, y: Array, w: Array) -> Array:
    return jac(t, y)@w

def step(...):
    ....
    y1 = ....
    y_error = ...
    w = update_w(t, y1, w)
    return ((y1, w), y_error, dense, ...)

I'm currently trying to wrap an arbitrary solver with some success like so:

def wrap_solver(solver_class):

    class WrappedSolver(solver_class):

        def update_aux(self, t0, x, aux):
             ....
             return aux

        def step(self, terms, t0, t1, y0, args, solver_state, made_jump):
            (y, y_error, dense_info, solver_state, solver_result) = super().step(terms, t0, t1, y0, args, solver_state, made_jump)
            x, aux = y
            new_aux = self.update_aux(t0, x, aux)
            return ((x, new_aux), y_error, dense_info, solver_state, solver_result)
      
    return WrappedSolver

But this seems very janky and I'm wondering if there's a cleaner way without exposing the step method. It's as if I need a custom term where one term is an ODETerm with control dt and the other is like a "StateTerm" also with control dt but it doesn't get integrated.

Start by reading the solver step API and the existing ODETerm handling described in the issue. Compare the wrapped-solver approach with the proposed control-dt state term that is updated but not integrated. Done means arbitrary solvers can maintain and update an auxiliary state without exposing or overriding step.