'Finds numerical value of a holonomic function using numerical methods.
(RK4 by default). A set of points (real or complex) must be provided
which will be the path for the numerical integration.
The path should be given as a list :math:`[x_1, x_2, ... x_n]`. The numerical
values will be computed at each point in this order
:math:`x_1 --> x_2 --> x_3 ... --> x_n`.
Returns values of the function at :math:`x_1, x_2, ... x_n` in a list.
Examples
>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols(\'x\')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),\'Dx\')
A straight line on the real axis from (0 to 1)
>>> r = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
Runge-Kutta 4th order on e^x from 0.1 to 1.
Exact solution at 1 is 2.71828182845905
>>> HolonomicFunction(Dx - 1, x, 0, [1]).evalf(r)
[1.10517083333333, 1.22140257085069, 1.34985849706254, 1.49182424008069,
1.64872063859684, 1.82211796209193, 2.01375162659678, 2.22553956329232,
2.45960141378007, 2.71827974413517]
Euler\'s method for the same
>>> HolonomicFunction(Dx - 1, x, 0, [1]).evalf(r, method=\'Euler\')
[1.1, 1.21, 1.331, 1.4641, 1.61051, 1.771561, 1.9487171, 2.14358881,
2.357947691, 2.5937424601]
One can also observe that the value obtained using Runge-Kutta 4th order
is much more accurate than Euler\'s method.'

def evalf(self, points, method='RK4', h=0.05, derivatives=False):

    from sympy.holonomic.numerical import _evalf
    lp = False
    if (not hasattr(points, '__iter__')):
        lp = True
        b = S(points)
        if (self.x0 == b):
            return _evalf(self, [b], method=method, derivatives=derivatives)[(-1)]
        if (not b.is_Number):
            raise NotImplementedError
        a = self.x0
        if (a > b):
            h = (- h)
        n = int(((b - a) / h))
        points = [(a + h)]
        for i in range((n - 1)):
            points.append((points[(-1)] + h))
    for i in roots(self.annihilator.parent.base.to_sympy(self.annihilator.listofpoly[(-1)]), self.x):
        if ((i == self.x0) or (i in points)):
            raise SingularityError(self, i)
    if lp:
        return _evalf(self, points, method=method, derivatives=derivatives)[(-1)]
    return _evalf(self, points, method=method, derivatives=derivatives)
