'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.'

def compute_residual(self, traj, ti, k=3, ext=2):

    soln = self.interpolate(traj, ti, k, 0, ext)
    deriv = self.interpolate(traj, ti, k, 1, ext)
    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))
    residual = (deriv - rhs_ode)
    return residual
