Function rkmk_dae_step_second_order_minimize 

pub fn rkmk_dae_step_second_order_minimize<Field, Tab, F, J, H, Z, T>(
    rate: &mut impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<Derivative<Field::Increment, T>, String>,
    function: impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<F, String>,
    jacobian: impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<J, String>,
    hessian: impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<H, String>,
    solver: &impl SecondOrderOptimization<F, J, H, Z>,
    point: &Field::Point,
    z: &Z,
    t: Quantity<T>,
    dt: Quantity<T>,
    scratch: &mut Vec<Field::Increment>,
    first_rate: Option<&Derivative<Field::Increment, T>>,
    equality_constraint: impl FnMut(Quantity<T>) -> EqualityConstraint,
    sparse: Option<SparseSolver>,
) -> Result<(Field::Point, Z, Option<Derivative<Field::Increment, T>>), IntegrationError>
where
    Field: Integrable,
    Tab: ButcherTableau,
    Field::Point: Clone,
    Field::Increment: Clone + Differentiable<T>,
    Z: Clone,
    T: Copy,
    Quantity<T>: Mul<Scalar, Output = Quantity<T>>,
    for<'a> &'a Derivative<Field::Increment, T>: Mul<Quantity<T>, Output = Field::Increment>,

rkmk_dae_step with the algebraic unknown resolved by second-order minimization at every stage abscissa, built from function/jacobian/ hessian/solver the same way rkmk_dae_step_first_order_root builds it for root-finding — the two are siblings so a model whose equilibrium is naturally posed as a potential (rather than a residual) gets the same manifold-aware stage-equilibrium step.