Function implicit_jacobian 

pub fn implicit_jacobian<F: Float>(
    tape: &mut BytecodeTape<F>,
    z_star: &[F],
    x: &[F],
    num_states: usize,
) -> Option<Vec<Vec<F>>>

Compute the full implicit Jacobian dz*/dx (m × n matrix).

Given a multi-output residual tape F: R^(m+n) → R^m with F(z*, x) = 0, computes dz*/dx = -F_z^{-1} · F_x via the Implicit Function Theorem.

The first num_states tape inputs are state variables z, the remaining are parameters x.

Returns None if F_z is singular.