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