Function solve_initial_condition_sensitivity 

pub fn solve_initial_condition_sensitivity<Sol, S, Sys>(
    system: &Sys,
    t0: S,
    tf: S,
    y0: &[S],
    options: &SolverOptions<S>,
) -> Result<StateTransitionResult<S>, SolverError>
where
    Sol: Solver<S>,
    S: Scalar,
    Sys: OdeSystem<S> + ?Sized,

Compute the state-transition matrix Φ(t) = ∂y(t) / ∂y₀ of the linearised flow of dy/dt = f(t, y), using any Solver.

Integrates the variational equations Ṡ = J_y · S with seed S(t₀) = I_N over the same AugmentedSystem used for parameter sensitivity, then extracts the state-transition trajectory.

The user supplies only the ODE system; the IC seeding (identity initial sensitivity, zero parameter forcing) is performed internally. J_y is obtained from OdeSystem::jacobian — the user’s analytical override if any, the trait’s forward-FD default otherwise. For an exact-to-round-off J_y from autodiff, wrap the system in [crate::AutodiffJacobianSystem] (behind the autodiff feature).

Example: scalar exponential decay

dy/dt = -k y has closed-form Φ(t) = exp(-k t).

use numra_ode::sensitivity::solve_initial_condition_sensitivity_with;
use numra_ode::{DoPri5, SolverOptions};

let k = 0.5_f64;
let result = solve_initial_condition_sensitivity_with::<DoPri5, f64, _>(
    move |_t, y, dy| { dy[0] = -k * y[0]; },
    &[1.0],
    0.0, 2.0,
    &SolverOptions::default().rtol(1e-9).atol(1e-12),
).unwrap();

let last = result.len() - 1;
let phi = result.phi_ij(last, 0, 0);
let analytical = (-k * result.t()[last]).exp();
assert!((phi - analytical).abs() < 1e-7);