Function implicit_euler_step 

pub fn implicit_euler_step(
    y: &[f64],
    jacobian: &[Vec<f64>],
    rhs: &[f64],
    dt: f64,
) -> Vec<f64>

Implicit Euler step: solve (I - dt*J) * y_new = y + dt*rhs.

Uses Gaussian elimination with partial pivoting.