pub fn solve_linear_system(a: &Matrix, b: &[f64]) -> Option<Vec<f64>>
Solve Ax = b via Gaussian elimination with partial pivoting.