Function thomas_solve

pub fn thomas_solve(
    lower: &[f64],
    diag: &[f64],
    upper: &[f64],
    rhs: &mut [f64],
) -> Vec<f64>

Solve a tridiagonal system (allocating version for tests/convenience).