docs: def compute_residual(self, traj, ti, k=3, ext=2):
          r"""
          The residual is the difference between the derivative of the B-spline
          approximation of the solution trajectory and the right-hand side of the
          original ODE evaluated along the approximated solution trajectory.
  
          Parameters
          ----------
          traj : array_like (float)
              Solution trajectory providing the data points for constructing the
              B-spline representation.
          ti : array_like (float)
              Array of values for the independent variable at which to
              interpolate the value of the B-spline.
          k : int, optional(default=3)
              Degree of the desired B-spline. Degree must satisfy
              :math:`1 \le k \le 5`.
          ext : int, optional(default=2)
              Controls the value of returned elements for outside the
              original knot sequence provided by traj. For extrapolation, set
              `ext=0`; `ext=1` returns zero; `ext=2` raises a `ValueError`.
  
          Returns
          -------
          residual : array (float)
              Difference between the derivative of the B-spline approximation
              of the solution trajectory and the right-hand side of the ODE
              evaluated along the approximated solution trajectory.
  
          """
          # B-spline approximations of the solution and its derivative
          soln = self.interpolate(traj, ti, k, 0, ext)
          deriv = self.interpolate(traj, ti, k, 1, ext)
  
          # rhs of ode evaluated along approximate solution
          T = ti.size
          rhs_ode = np.vstack(self.f(ti[i], soln[i, 1:], *self.f_params)
                              for i in range(T))
          rhs_ode = np.hstack((ti[:, np.newaxis], rhs_ode))
  
          # should be roughly zero everywhere (if approximation is any good!)
          residual = deriv - rhs_ode
  
          return residual