code: def collocation_fun(fun, y, p, x, h):
      """Evaluate collocation residuals.
  
      This function lies in the core of the method. The solution is sought
      as a cubic C1 continuous spline with derivatives matching the ODE rhs
      at given nodes `x`. Collocation conditions are formed from the equality
      of the spline derivatives and rhs of the ODE system in the middle points
      between nodes.
  
      Such method is classified to Lobbato IIIA family in ODE literature.
      Refer to [1]_ for the formula and some discussion.
  
      Returns
      -------
      col_res : ndarray, shape (n, m - 1)
          Collocation residuals at the middle points of the mesh intervals.
      y_middle : ndarray, shape (n, m - 1)
          Values of the cubic spline evaluated at the middle points of the mesh
          intervals.
      f : ndarray, shape (n, m)
          RHS of the ODE system evaluated at the mesh nodes.
      f_middle : ndarray, shape (n, m - 1)
          RHS of the ODE system evaluated at the middle points of the mesh
          intervals (and using `y_middle`).
  
      References
      ----------
      .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
             Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
             Number 3, pp. 299-316, 2001.
      """
      f = fun(x, y, p)
      y_middle = 0.5 * (y[:, 1:] + y[:, :-1]) - 0.125 * h * (
          f[:, 1:] - f[:, :-1]
      )
      f_middle = fun(x[:-1] + 0.5 * h, y_middle, p)
      col_res = (
          y[:, 1:]
          - y[:, :-1]
          - h / 6 * (f[:, :-1] + f[:, 1:] + 4 * f_middle)
      )
  
      return col_res, y_middle, f, f_middle
  
imports: from __future__ import division, print_function, absolute_import
  from warnings import warn
  import numpy
  from numpy.linalg import norm, pinv
  from scipy.sparse import coo_matrix, csc_matrix
  from scipy.sparse.linalg import splu
  from scipy.optimize import OptimizeResult
label_desc: Use image node_id 6 for calling a global function with example usage: collocation_fun(fun, y, p, x, h) and returns: col_res,  y_middle,  f,  f_middle