// ode_solver.mbt — ODE solvers for Neural ODEs (v0.90.0).
//
// Integrators that advance the state x(t) of an ODE dx/dt = f(x, t; θ)
// over a time interval [t0, t1] using a fixed step size `h` and a
// fixed number of `n_steps`. Both solvers return:
//   - the final state x(t1)
//   - the full state trajectory flat `[n_steps+1 × state_dim]` for use
//     by the adjoint method (v0.91.0).
//
// Scope of v0.90.0:
//   - ode_euler_solve: explicit Euler integration (1st order)
//   - ode_heun_solve: Heun / improved Euler / RK2 (2nd order)
//
// Reference: Hairer & Wanner 1991 "Solving Ordinary Differential
// Equations"; Chen et al. 2018 (Neural ODEs) uses these solvers.

///|
/// Explicit Euler integration. Step: x_{n+1} = x_n + h · f(x_n, t_n).
/// O(h) local error, O(h) global error. Cheapest per step but
/// accumulates error quickly for stiff systems.
pub fn ode_euler_solve(
  func : ODEFunc,
  x0 : Array[Float],
  t0 : Float,
  t1 : Float,
  n_steps : Int,
) -> (Array[Float], Array[Float]) {
  let state_dim = func.state_dim
  let trajectory : Array[Float] = Array::make(
    (n_steps + 1) * state_dim, 0.0F,
  )
  // Initial state.
  let x : Array[Float] = Array::make(state_dim, 0.0F)
  for i in 0.. (Array[Float], Array[Float]) {
  let state_dim = func.state_dim
  let trajectory : Array[Float] = Array::make(
    (n_steps + 1) * state_dim, 0.0F,
  )
  let x : Array[Float] = Array::make(state_dim, 0.0F)
  for i in 0..