pub fn controllability_matrix(a: &[Vec<f64>], b: &[Vec<f64>]) -> Vec<Vec<f64>>
Compute controllability matrix: [B, AB, A²B, …, A^(n-1)B].