Iterate over matrices of varying shape

Question

I'm implementing MATLAB's ode45 in Rust, but attempting to make a major improvement. If a user of MATLAB's ode45 wants to integrate a matrix and a vector simultaneously, they have to flatten the matrix and append or prepend the vector before passing it to MATLAB's ode45 only to undo that process when attempting to make use of the solution. I would like for a user of my ode45 to be able to pass a collection¹ containing arbitrarily-shaped nalgebra Matrix types.

I've seen that I can't iterate over tuples unless I use a fancy macro and implement a trait for every size of tuple the macro is to be used on. I also can't use a Matrix array, because the types aren't strictly the same if the shapes are arbitrary.

An important detail about integrators is that you must provide a function that models the first-order differential equation. From the user's perspective, this would be fairly easy, since they could simply write a function that takes their specific tuple and returns the same. The trouble is figuring out how to handle whatever the user may want to integrate simultaneously.

Below is a snippet that illustrates what I want to be able to do:

let system_dynamics(y: &(&SVector<f64, 4>, &SMatrix<f64, 4, 4>)
    -> (SVector<f64, 4>, SMatrix<f64, 4, 4>) {
  return (SVector::<f64, 4>::new(y.0[2], y.0[3], 0.0, -9.81),
          SMatrix::<f64, 4, 4>::new(y.0[2],    0.0,    0.0,    0.0,
                                       0.0, y.0[2],    0.0,    0.0,
                                       0.0,    0.0, y.0[2],    0.0,
                                       0.0,    0.0,    0.0, y.0[2])
  );
}

let y0 = (SVector::<f64, 4>::new(0.0, 0.0, 10.0, 10.0),
          SMatrix::<f64, 4, 4>::zeros());
let t0 = 0.0;
let tf = 10.0;
let (t, y) = ode45(|_t, y| system_dynamics(&y), [t0, tf], y0);

println!("{}: ({}, {})", t.last().unwrap, y.last().unwrap().0, y.last().unwrap().1);
// Should print something along the lines of
// 10: (
//   ┌        ┐
//   │    100 │
//   │ -390.5 │
//   │     10 │
//   │  -88.1 │
//   └        ┘
// ,
//   ┌                    ┐
//   │ 100  0.0  0.0  0.0 │
//   │ 0.0  100  0.0  0.0 │
//   │ 0.0  0.0  100  0.0 │
//   │ 0.0  0.0  0.0  100 │
//   └                    ┘
// )

I already have a working and tested implementation of ode45 that works on an arbitrarily-shaped Matrix; all I need help with is figuring out how to accept multiple Matrix types of arbitrary dimension that need to be integrated simultaneously.

---

^{1: by "collection", I mean something of arbitrary length that can be iterated over.}