docs: def solution(self, x0, y0):
          """ Create a solution function ``y(x)`` such that ``y(x0) = y0``.
  
          A list of solution values ``[y(x0), y(x1) ...]`` is returned if the
          function is called with a list ``[x0, x1 ...]`` of ``x`` values.
          """
          def soln(x):
              if numpy.size(x) > 1:
                  x = [soln.x] + list(x)
                  ans = self(soln.y, interval=x)
                  soln.x = x[-1]
                  soln.y = ans[-1]
                  return ans
              else:
                  soln.y = self(soln.y, interval=(soln.x, x))
                  soln.x = x
                  return soln.y
          soln.x = x0
          soln.y = y0
          return soln