pub fn solve_linear_system( a: &Array2<f64>, b: &Array1<f64>, ) -> IntegrateResult<Array1<f64>>
Solve linear system Ax = b using LU decomposition with partial pivoting