///|
type OdeSystemFn = (@mpf.RawMpf, ArrayView[@mpf.RawMpf]) -> Array[@mpf.RawMpf]

///|
fn MPContext::p3_ode_expect_dim(
  self : MPContext,
  name : String,
  y : ArrayView[@mpf.RawMpf],
  n : Int,
) -> Unit raise MPError {
  ignore(self)
  if y.length() != n {
    raise ValueError("\{name}: derivative dimension mismatch")
  }
}

///|
fn MPContext::p3_ode_shift(
  self : MPContext,
  y : ArrayView[@mpf.RawMpf],
  k : ArrayView[@mpf.RawMpf],
  scale : @mpf.RawMpf,
  prec : Int,
) -> Array[@mpf.RawMpf] {
  let out : Array[@mpf.RawMpf] = []
  for i in 0.. Array[@mpf.RawMpf] raise MPError {
  let h6 = p2_mpf_div_int(h, 6, prec, @mpf.round_nearest)
  let out : Array[@mpf.RawMpf] = []
  for i in 0.. Array[@mpf.RawMpf] raise MPError {
  let n = y.length()
  let p = self.p2_work_prec()
  let h2 = p2_mpf_div_int(h, 2, p, @mpf.round_nearest)
  let t2 = @mpf.mpf_add(t, h2, p, @mpf.round_nearest)
  let t4 = @mpf.mpf_add(t, h, p, @mpf.round_nearest)

  let k1 = f(t, y)
  self.p3_ode_expect_dim("ode_step_rk4", k1, n)
  let y2 = self.p3_ode_shift(y, k1, h2, p)

  let k2 = f(t2, y2)
  self.p3_ode_expect_dim("ode_step_rk4", k2, n)
  let y3 = self.p3_ode_shift(y, k2, h2, p)

  let k3 = f(t2, y3)
  self.p3_ode_expect_dim("ode_step_rk4", k3, n)
  let y4 = self.p3_ode_shift(y, k3, h, p)

  let k4 = f(t4, y4)
  self.p3_ode_expect_dim("ode_step_rk4", k4, n)

  self.p3_ode_finalize(y, k1, k2, k3, k4, h, p)
}

///|
pub fn MPContext::ode_solve(
  self : MPContext,
  f : OdeSystemFn,
  t0 : @mpf.RawMpf,
  y0 : ArrayView[@mpf.RawMpf],
  t1 : @mpf.RawMpf,
  steps? : Int = 200,
) -> Array[@mpf.RawMpf] raise MPError {
  if steps < 1 {
    raise ValueError("ode_solve: steps must be positive")
  }
  let p = self.p2_work_prec()
  let h = p2_mpf_div(
    @mpf.mpf_sub(t1, t0, p, @mpf.round_nearest),
    @mpf.from_int(steps),
    p,
    @mpf.round_nearest,
  )
  let mut t = t0
  let mut y = p3_copy_vec(y0)
  for _ in 0.. Array[Array[@mpf.RawMpf]] raise MPError {
  if steps_per_interval < 1 {
    raise ValueError("ode_solve_grid: steps_per_interval must be positive")
  }
  if ts.length() == 0 {
    return []
  }
  if !@mpf.mpf_eq(ts[0], t0) {
    raise ValueError("ode_solve_grid: ts[0] must equal t0")
  }
  let out : Array[Array[@mpf.RawMpf]] = []
  let mut t = t0
  let mut y = p3_copy_vec(y0)
  out.push(p3_copy_vec(y))
  for i in 1.. (Array[@mpf.RawMpf], MpfMatrix) raise MPError {
  if a.rows != a.cols {
    raise ValueError("eigsy: matrix must be square")
  }
  let n = a.rows
  if n == 0 {
    return ([], { rows: 0, cols: 0, data: [] })
  }
  if n == 1 {
    return (
      [@mpf.mpf_pos(a.data[0], self.precision(), self.round_mode())],
      self.eye(1),
    )
  }
  for i in 0.. @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => {
      let s = if self.precision() > 0 { -(self.precision() / 2) } else { -40 }
      @mpf.from_man_exp(1N, s, 0, @mpf.round_down)
    }
  }
  let steps = if max_steps > 0 { max_steps } else { n * n * 32 }
  let work = p3_copy_vec(a.data)
  let q0 = self.eye(n)
  let q = p3_copy_vec(q0.data)
  let mut done = false
  for _ in 0.. raise from_mpf_error(err)
    }
    let denom = @mpf.mpf_add(tau_abs, root, p, @mpf.round_nearest)
    let mut t = if @mpf.is_zero(denom) {
      @mpf.fone
    } else {
      p2_mpf_div(@mpf.fone, denom, p, @mpf.round_nearest)
    }
    if tau.sign == 1 {
      t = @mpf.mpf_neg(t, p, @mpf.round_nearest)
    }
    let c = p2_mpf_div(
      @mpf.fone,
      @mpf.mpf_sqrt(
        @mpf.mpf_add(
          @mpf.fone,
          @mpf.mpf_mul(t, t, p, @mpf.round_nearest),
          p,
          @mpf.round_nearest,
        ),
        p,
        @mpf.round_nearest,
      ) catch {
        err => raise from_mpf_error(err)
      },
      p,
      @mpf.round_nearest,
    )
    let s = @mpf.mpf_mul(t, c, p, @mpf.round_nearest)
    let c2 = @mpf.mpf_mul(c, c, p, @mpf.round_nearest)
    let s2 = @mpf.mpf_mul(s, s, p, @mpf.round_nearest)
    let sc2 = @mpf.mpf_mul_int(
      @mpf.mpf_mul(s, c, p, @mpf.round_nearest),
      2,
      p,
      @mpf.round_nearest,
    )
    for k in 0..