Function compute_consistent_initial_tol 

pub fn compute_consistent_initial_tol<S, Sys>(
    system: &Sys,
    t0: S,
    y0: &[S],
    tol: S,
    max_iter: usize,
) -> Result<Vec<S>, String>
where
    S: Scalar + SimpleEntity + Conjugate<Canonical = S> + ComplexField,
    Sys: OdeSystem<S>,

Compute consistent initial conditions with user-specified tolerance.

Algorithm

Uses Newton iteration with a full Jacobian on the algebraic variable subset:

1. Keep differential variables fixed
2. Evaluate residual g(t₀, y) on algebraic equations
3. Build full Jacobian ∂g/∂y_alg via finite differences
4. Solve J_alg * Δy_alg = -g using LU factorization
5. Update y_alg += Δy_alg
6. Repeat until convergence

This handles coupled algebraic constraints correctly (unlike a diagonal Jacobian approach which fails when constraints share variables).