Add RK4 solver for OdeProblem

[dyreby]

Forward Euler landed in #242. This issue tracks the RK4 integrator.

Add transient::rk4 following the same structure as transient::euler:

• solve / solve_unobserved API
• Own Config, Solution, Error, Event, Action types
• Event emitted after each accepted step
• Action::StopEarly for observer-controlled termination

Design note

RK4 uses four intermediate evaluations per step (k1–k4), none of which should trigger finalize_step — that hook runs only on the accepted step, same as Euler.

References

• Add transient stepping solvers for OdeProblem #237 (original, also covered Euler)
• Add forward Euler solver for OdeProblem #242 (Euler implementation, for reference)

[dyreby]

Design notes from pre-implementation review

RK4 requires more from Delta and Derivative types than Euler does, because of fractional deltas and derivative arithmetic.

What RK4 needs that Euler doesn't

k1 = f(t,        y_n)
k2 = f(t + dt/2, y_n + k1 * dt/2)
k3 = f(t + dt/2, y_n + k2 * dt/2)
k4 = f(t + dt,   y_n + k3 * dt)

y_{n+1} = y_n + (k1 + 2*k2 + 2*k3 + k4) * dt/6

1. Fractional deltas — dt/2 and dt/6 needed for intermediate build_input calls. Delta has no arithmetic in OdeProblem or StepIntegrable as currently defined.

2. Derivative arithmetic — combining k1 + 2*k2 + 2*k3 + k4 requires Add and Mul<f64> on the derivative type. StepIntegrable::step takes one derivative at a time and says nothing about arithmetic between derivatives.

3. Intermediate build_input calls — for k2, build_input(base, intermediate_state, dt/2) must receive dt/2 not dt, or time tracking in the input will be wrong.

Approach: explicit bounds on solve, with alias

Keep OdeProblem and StepIntegrable as-is. RK4's solve advertises what it needs via trait bounds. To keep the call site clean, export aliases with blanket impls:

pub trait Rk4Delta: Copy + Mul<f64, Output = Self> {}
impl<T: Copy + Mul<f64, Output = T>> Rk4Delta for T {}

pub trait Rk4Derivative: Add<Output = Self> + Mul<f64, Output = Self> + Clone {}
impl<T: Add<Output = T> + Mul<f64, Output = T> + Clone> Rk4Derivative for T {}

Blanket impls mean nothing extra to implement — any type satisfying the bounds gets the alias for free. solve then reads:

pub fn solve<M, P, Obs>(...)
where
    M: Model,
    P: OdeProblem<Input = M::Input, Output = M::Output>,
    P::Delta: Rk4Delta,
    DerivativeOf<P::State, P::Delta>: Rk4Derivative,
    ...

This is honest about requirements, clean at the call site, and leaves room to tighten or loosen bounds later without a trait redesign.