Repository: sympy/sympy

Issue:
Add a docstring how to plot a holonomic function
Something like this:
```
import sympy.holonomic
import numpy as np
r = np.linspace(1, 5, 100)
y = sympy.holonomic.from_sympy(sin(x)**2/x).evalf(r)
plot(r, y)
```

This should be part of some larger docstring, or part of a tutorial for the holonomic module. Here is a jupyter notebook showing that it works:

http://nbviewer.jupyter.org/gist/certik/849a62d23f615d89dd37585b76159e78

cc @shubhamtibra .

Discussion:

Patch:
diff --git a/setup.py b/setup.py
index 097efc5f42..bee2324d13 100755
--- a/setup.py
+++ b/setup.py
@@ -86,6 +86,7 @@
     'sympy.functions.special',
     'sympy.functions.special.benchmarks',
     'sympy.geometry',
+    'sympy.holonomic',
     'sympy.integrals',
     'sympy.integrals.benchmarks',
     'sympy.interactive',
@@ -270,6 +271,7 @@ def run(self):
     'sympy.functions.elementary.tests',
     'sympy.functions.special.tests',
     'sympy.geometry.tests',
+    'sympy.holonomic.tests',
     'sympy.integrals.tests',
     'sympy.interactive.tests',
     'sympy.liealgebras.tests',
diff --git a/sympy/holonomic/holonomic.py b/sympy/holonomic/holonomic.py
index 2a2c1ef65e..3868cb4e80 100644
--- a/sympy/holonomic/holonomic.py
+++ b/sympy/holonomic/holonomic.py
@@ -2,7 +2,8 @@
 
 from __future__ import print_function, division
 
-from sympy import symbols, Symbol, diff, S, Dummy, Order, rf, meijerint, I, solve
+from sympy import (symbols, Symbol, diff, S, Dummy, Order, rf, meijerint, I,
+    solve, limit)
 from sympy.printing import sstr
 from .linearsolver import NewMatrix
 from .recurrence import HolonomicSequence, RecurrenceOperator, RecurrenceOperators
@@ -323,6 +324,14 @@ def __eq__(self, other):
             else:
                 return False
 
+    def is_singular(self, x0=0):
+        """
+        Checks if the differential equation is singular at x0.
+        """
+
+        domain = str(self.parent.base.domain)[0]
+        return x0 in roots(self.listofpoly[-1].rep, filter=domain)
+
 
 class HolonomicFunction(object):
     """
@@ -336,6 +345,21 @@ class HolonomicFunction(object):
     Given two Holonomic Functions f and g, their sum, product,
     integral and derivative is also a Holonomic Function.
 
+    To plot a Holonomic Function, one can use `.evalf()` for numerical computation.
+    Here's an example on `sin(x)**2/x` using numpy and matplotlib.
+
+    ``
+    import sympy.holonomic
+    from sympy import var, sin
+    import matplotlib.pyplot as plt
+    import numpy as np
+    var("x")
+    r = np.linspace(1, 5, 100)
+    y = sympy.holonomic.from_sympy(sin(x)**2/x, x0=1).evalf(r)
+    plt.plot(r, y, label="holonomic function")
+    plt.show()
+    ``
+
     Examples
     ========
 
@@ -482,6 +506,7 @@ def __add__(self, other):
         sol = _normalize(sol.listofpoly, self.annihilator.parent, negative=False)
         # solving initial conditions
         if self._have_init_cond and other._have_init_cond:
+
             if self.x0 == other.x0:
                 # try to extended the initial conditions
                 # using the annihilator
@@ -489,14 +514,29 @@ def __add__(self, other):
                 y0_other = _extend_y0(other, sol.order)
                 y0 = [a + b for a, b in zip(y0_self, y0_other)]
                 return HolonomicFunction(sol, self.x, self.x0, y0)
+
             else:
-                # initial conditions for different points
-                # to be implemented
-                pass
+
+                if self.x0 == 0 and not self.annihilator.is_singular() \
+                and not other.annihilator.is_singular():
+                    return self + other.change_ics(0)
+
+                elif other.x0 == 0 and not self.annihilator.is_singular() \
+                and not other.annihilator.is_singular():
+                    return self.change_ics(0) + other
+
+                else:
+
+                    if not self.annihilator.is_singular(x0=self.x0) and \
+                    not other.annihilator.is_singular(x0=self.x0):
+                        return self + other.change_ics(self.x0)
+
+                    else:
+                        return self.change_ics(other.x0) + other
 
         return HolonomicFunction(sol, self.x)
 
-    def integrate(self, *args):
+    def integrate(self, limits):
         """
         Integrate the given holonomic function. Limits can be provided,
         Initial conditions can only be computed when limits are (x0, x).
@@ -518,19 +558,126 @@ def integrate(self, *args):
         HolonomicFunction((1)Dx + (1)Dx**3, x), f(0) = 0, f'(0) = 1, f''(0) = 0
         """
 
-        # just multiply by Dx from right
+        # to get the annihilator, just multiply by Dx from right
         D = self.annihilator.parent.derivative_operator
-        if (not args) or (not self._have_init_cond):
+
+        # for indefinite integration
+        if (not limits) or (not self._have_init_cond):
             return HolonomicFunction(self.annihilator * D, self.x)
 
-        # definite integral if limits are (x0, x)
-        if len(args) == 1 and len(args[0]) == 3 and args[0][1] == self.x0 and args[0][2] == self.x:
+        # definite integral
+        # initial conditions for the answer will be stored at point `a`,
+        # where `a` is the lower limit of the integrand
+        if hasattr(limits, "__iter__"):
+
+            if len(limits) == 3 and limits[0] == self.x:
+                x0 = self.x0
+                a = limits[1]
+                b = limits[2]
+
+        else:
+            x0 = self.x0
+            a = self.x0
+            b = self.x
+
+        if x0 == a:
             y0 = [S(0)]
             y0 += self.y0
-            return HolonomicFunction(self.annihilator * D, self.x, self.x0, y0)
+
+        # use evalf to get the values at `a`
+        else:
+            y0 = [S(0)]
+            tempy0 = self.evalf(a, derivatives=True)
+            y0 += tempy0
+
+        # if the upper limit is `x`, the answer will be a function
+        if b == self.x:
+            return HolonomicFunction(self.annihilator * D, self.x, a, y0)
+
+        # if the upper limits is a Number, a numerical value will be returned
+        elif S(b).is_Number:
+            return HolonomicFunction(self.annihilator * D, self.x, a, y0).evalf(b)
 
         return HolonomicFunction(self.annihilator * D, self.x)
 
+    def diff(self, *args):
+        """
+        Differentiation of the given Holonomic function.
+
+        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(ZZ.old_poly_ring(x),'Dx')
+
+        # derivative of sin(x)
+        >>> HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]).diff().to_sympy()
+        cos(x)
+
+        # derivative of e^2*x
+        >>> HolonomicFunction(Dx - 2, x, 0, 1).diff().to_sympy()
+        2*exp(2*x)
+
+        See Also
+        =======
+
+        .integrate()
+        """
+
+        if args:
+            if args[0] != self.x:
+                return S(0)
+            elif len(args) == 2:
+                sol = self
+                for i in range(args[1]):
+                    sol = sol.diff(args[0])
+                return sol
+
+        ann = self.annihilator
+        dx = ann.parent.derivative_operator
+
+        # if the function is constant.
+        if ann.listofpoly[0] == ann.parent.base.zero and ann.order == 1:
+            return S(0)
+
+        # if the coefficient of y in the differential equation is zero.
+        # a shifting is done to compute the answer in this case.
+        elif ann.listofpoly[0] == ann.parent.base.zero:
+
+            sol = DifferentialOperator(ann.listofpoly[1:], ann.parent)
+
+            if self._have_init_cond:
+                return HolonomicFunction(sol, self.x, self.x0, self.y0[1:])
+
+            else:
+                return HolonomicFunction(sol, self.x)
+
+        # the general algorithm
+        R = ann.parent.base
+        K = R.get_field()
+
+        seq_dmf = [K.new(i.rep) for i in ann.listofpoly]
+
+        # -y = a1*y'/a0 + a2*y''/a0 ... + an*y^n/a0
+        rhs = [i / seq_dmf[0] for i in seq_dmf[1:]]
+        rhs.insert(0, K.zero)
+
+        # differentiate both lhs and rhs
+        sol = _derivate_diff_eq(rhs)
+
+        # add the term y' in lhs to rhs
+        sol = _add_lists(sol, [K.zero, K.one])
+
+        sol = _normalize(sol[1:], self.annihilator.parent, negative=False)
+
+        if not self._have_init_cond:
+            return HolonomicFunction(sol, self.x)
+
+        y0 = _extend_y0(self, sol.order + 1)[1:]
+        return HolonomicFunction(sol, self.x, self.x0, y0)
 
     def __eq__(self, other):
         if self.annihilator == other.annihilator:
@@ -635,6 +782,8 @@ def __mul__(self, other):
         sol_ann = _normalize(sol[0][0:], self.annihilator.parent, negative=False)
 
         if self._have_init_cond and other._have_init_cond:
+
+            # if both the conditions are at same point
             if self.x0 == other.x0:
 
                 # try to find more inital conditions
@@ -657,9 +806,28 @@ def __mul__(self, other):
                             sol += coeff[j][k]* y0_self[j] * y0_other[k]
 
                     y0.append(sol)
+
                 return HolonomicFunction(sol_ann, self.x, self.x0, y0)
+
+            # if the points are different, consider one
             else:
-                raise NotImplementedError
+
+                if self.x0 == 0 and not ann_self.is_singular() \
+                and not ann_other.is_singular():
+                    return self * other.change_ics(0)
+
+                elif other.x0 == 0 and not ann_self.is_singular() \
+                and not ann_other.is_singular():
+                    return self.change_ics(0) * other
+
+                else:
+
+                    if not ann_self.is_singular(x0=self.x0) and \
+                    not ann_other.is_singular(x0=self.x0):
+                        return self * other.change_ics(self.x0)
+
+                    else:
+                        return self.change_ics(other.x0) * other
 
         return HolonomicFunction(sol_ann, self.x)
 
@@ -770,10 +938,15 @@ def composition(self, expr, *args):
     def to_sequence(self, lb=True):
         """
         Finds the recurrence relation in power series expansion
-        of the function.
+        of the function about origin.
+
+        Returns a tuple (R, n0), with R being the Recurrence relation
+        and a non-negative integer `n0` such that the recurrence relation
+        holds for all n >= n0.
 
-        Returns a tuple of the recurrence relation and a value `n0` such that
-        the recurrence relation holds for all n >= n0.
+        If it's not possible to numerically compute a initial condition,
+        it is returned as a symbol C_j, denoting the coefficient of x^j
+        in the power series about origin.
 
         Examples
         ========
@@ -858,12 +1031,6 @@ def to_sequence(self, lb=True):
         # the recurrence relation
         sol = RecurrenceOperator(sol, R)
 
-        if not self._have_init_cond or self.x0 != 0:
-            if lb:
-                return (HolonomicSequence(sol), smallest_n)
-            else:
-                return HolonomicSequence(sol)
-
         # computing the initial conditions for recurrence
         order = sol.order
         all_roots = roots(sol.listofpoly[-1].rep, filter='Z')
@@ -884,34 +1051,61 @@ def to_sequence(self, lb=True):
 
         # if sufficient conditions can't be computed then
         # try to use the series method i.e.
-        # equate the coefficients of x^k in series to zero.
+        # equate the coefficients of x^k in the equation formed by
+        # substituting the series in differential equation, to zero.
         if len(u0) < order:
 
             for i in range(degree):
                 eq = S(0)
 
                 for j in dict1:
-                    if i + j[0] < len(u0):
+
+                    if i + j[0] < 0:
+                        dummys[i + j[0]] = S(0)
+
+                    elif i + j[0] < len(u0):
                         dummys[i + j[0]] = u0[i + j[0]]
 
                     elif not i + j[0] in dummys:
-                        dummys[i + j[0]] = Dummy('x_%s' %(i + j[0]))
+                        dummys[i + j[0]] = Symbol('C_%s' %(i + j[0]))
 
                     if j[1] <= i:
                         eq += dict1[j].subs(n, i) * dummys[i + j[0]]
 
                 eqs.append(eq)
 
+            # solve the system of equations formed
             soleqs = solve(eqs)
 
+            if isinstance(soleqs, dict):
+
+                for i in range(len(u0), order):
+
+                    if i not in dummys:
+                        dummys[i] = Symbol('C_%s' %i)
+
+                    if dummys[i] in soleqs:
+                        u0.append(soleqs[dummys[i]])
+
+                    else:
+                        u0.append(dummys[i])
+
+                if lb:
+                    return (HolonomicSequence(sol, u0), smallest_n)
+                return HolonomicSequence(sol, u0)
+
             for i in range(len(u0), order):
+
+                if i not in dummys:
+                    dummys[i] = Symbol('C_%s' %i)
+
                 s = False
                 for j in soleqs:
                     if dummys[i] in j:
                         u0.append(j[dummys[i]])
                         s = True
                 if not s:
-                    break
+                    u0.append(dummys[i])
 
         if lb:
             return (HolonomicSequence(sol, u0), smallest_n)
@@ -1022,7 +1216,7 @@ def _pole_degree(poly):
             y *= x - i
         return roots(s.rep, filter='R').keys()
 
-    def evalf(self, points, method='RK4'):
+    def evalf(self, points, method='RK4', h=0.05, derivatives=False):
         """
         Finds numerical value of a holonomic function using numerical methods.
         (RK4 by default). A set of points (real or complex) must be provided
@@ -1062,10 +1256,33 @@ def evalf(self, points, method='RK4'):
         """
 
         from sympy.holonomic.numerical import _evalf
+        lp = False
+
+        # if a point `b` is given instead of a mesh
+        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.listofpoly[-1].rep):
             if i == self.x0 or i in points:
                 raise TypeError("Provided path contains a singularity")
-        return _evalf(self, points, method=method)
+
+        if lp:
+            return _evalf(self, points, method=method, derivatives=derivatives)[-1]
+        return _evalf(self, points, method=method, derivatives=derivatives)
 
     def _subs(self, z):
         """
@@ -1123,6 +1340,18 @@ def to_hyper(self):
         # order of the recurrence relation
         m = r.order
 
+        # when no recurrence exists, and the power series have finite terms
+        if m == 0:
+            nonzeroterms = roots(r.listofpoly[0].rep, filter='Z')
+            if max(nonzeroterms) >= len(u0):
+                raise NotImplementedError("Initial conditions aren't sufficient")
+            sol = S(0)
+            for i in nonzeroterms:
+                if i >= 0:
+                    sol += u0[i] * x**i
+            return sol
+
+
         if smallest_n + m > len(u0):
             raise NotImplementedError("Can't compute sufficient Initial Conditions")
 
@@ -1135,7 +1364,7 @@ def to_hyper(self):
                 is_hyper = False
                 break
 
-        if not is_hyper or m == 0:
+        if not is_hyper:
             raise TypeError("The series is not Hypergeometric")
 
         a = r.listofpoly[0]
@@ -1214,6 +1443,37 @@ def to_sympy(self):
 
         return hyperexpand(self.to_hyper()).simplify()
 
+    def change_ics(self, b):
+        """
+        Changes the point `x0` to `b` for initial conditions.
+
+        Examples
+        ========
+
+        >>> from sympy.holonomic import from_sympy
+        >>> from sympy import symbols, sin, cos, exp
+        >>> x = symbols('x')
+
+        >>> from_sympy(sin(x)).change_ics(1)
+        HolonomicFunction((1) + (1)Dx**2, x), f(1) = sin(1), f'(1) = cos(1)
+
+        >>> from_sympy(exp(x)).change_ics(2)
+        HolonomicFunction((-1) + (1)Dx, x), f(2) = exp(2)
+        """
+
+        symbolic = True
+
+        try:
+            sol = from_sympy(self.to_sympy(), x0=b)
+        except:
+            symbolic = False
+
+        if symbolic:
+            return sol
+
+        y0 = self.evalf(b, derivatives=True)
+        return HolonomicFunction(self.annihilator, self.x, b, y0)
+
 
 def from_hyper(func, x0=0, evalf=False):
     """
@@ -1376,7 +1636,7 @@ def _find_conditions(simp, x, x0, order, evalf=False):
 from sympy.integrals.meijerint import _mytype
 
 
-def from_sympy(func, x=None, initcond=True):
+def from_sympy(func, x=None, initcond=True, x0=0):
     """
     Uses `meijerint._rewrite1` to convert to `meijerg` function and then
     eventually to Holonomic Functions. Only works when `meijerint._rewrite1`
@@ -1404,7 +1664,7 @@ def from_sympy(func, x=None, initcond=True):
         x = func.atoms(Symbol).pop()
 
     # try to convert if the function is polynomial or rational
-    solpoly = _convert_poly_rat(func, x, initcond=initcond)
+    solpoly = _convert_poly_rat(func, x, initcond=initcond, x0=x0)
     if solpoly:
         return solpoly
 
@@ -1425,7 +1685,6 @@ def from_sympy(func, x=None, initcond=True):
             sol = _convert_meijerint(func, x, initcond=False)
             if not sol:
                 raise NotImplementedError
-            x0 = 0
             y0 = _find_conditions(func, x, x0, sol.annihilator.order)
             while not y0:
                 x0 += 1
@@ -1435,7 +1694,6 @@ def from_sympy(func, x=None, initcond=True):
         if not initcond:
             return sol.composition(func.args[0])
 
-        x0 = 0
         y0 = _find_conditions(func, x, x0, sol.annihilator.order)
         while not y0:
             x0 += 1
@@ -1464,7 +1722,6 @@ def from_sympy(func, x=None, initcond=True):
     if not initcond:
         return sol
 
-    x0 = 0
     y0 = _find_conditions(func, x, x0, sol.annihilator.order)
     while not y0:
         x0 += 1
@@ -1484,8 +1741,7 @@ def _normalize(list_of, parent, negative=True):
     denom = []
     base = parent.base
     K = base.get_field()
-    R = ZZ.old_poly_ring(base.gens[0])
-    lcm_denom = R.from_sympy(S(1))
+    lcm_denom = base.from_sympy(S(1))
     list_of_coeff = []
 
     # convert polynomials to the elements of associated
@@ -1499,15 +1755,10 @@ def _normalize(list_of, parent, negative=True):
             list_of_coeff.append(j)
 
         # corresponding numerators of the sequence of polynomials
-        num.append(base(list_of_coeff[i].num))
+        num.append(list_of_coeff[i].numer())
 
         # corresponding denominators
-        den = list_of_coeff[i].den
-        if isinstance(den[0], PythonRational):
-            for i, j in enumerate(den):
-                den[i] = j.p
-
-        denom.append(R(den))
+        denom.append(list_of_coeff[i].denom())
 
     # lcm of denominators in the coefficients
     for i in denom:
@@ -1522,7 +1773,7 @@ def _normalize(list_of, parent, negative=True):
     for i, j in enumerate(list_of_coeff):
         list_of_coeff[i] = j * lcm_denom
 
-    gcd_numer = base.from_FractionField(list_of_coeff[-1], K)
+    gcd_numer = base((list_of_coeff[-1].numer() / list_of_coeff[-1].denom()).rep)
 
     # gcd of numerators in the coefficients
     for i in num:
@@ -1532,7 +1783,8 @@ def _normalize(list_of, parent, negative=True):
 
     # divide all the coefficients by the gcd
     for i, j in enumerate(list_of_coeff):
-        list_of_coeff[i] = base.from_FractionField(j / gcd_numer, K)
+        frac_ans = j / gcd_numer
+        list_of_coeff[i] = base((frac_ans.numer() / frac_ans.denom()).rep)
 
     return DifferentialOperator(list_of_coeff, parent)
 
@@ -1618,7 +1870,7 @@ def DMFdiff(frac):
     return K((sol_num.rep, sol_denom.rep))
 
 
-def DMFsubs(frac, x0):
+def DMFsubs(frac, x0, mpm=False):
     # substitute the point x0 in DMF object of the form p/q
     if not isinstance(frac, DMF):
         return frac
@@ -1628,20 +1880,23 @@ def DMFsubs(frac, x0):
     sol_p = S(0)
     sol_q = S(0)
 
+    if mpm:
+        from mpmath import mp
+
     for i, j in enumerate(reversed(p)):
-        if isinstance(j, PythonRational):
-            j = sympify(j)
+        if mpm:
+            j = sympify(j)._to_mpmath(mp.prec)
         sol_p += j * x0**i
 
     for i, j in enumerate(reversed(q)):
-        if isinstance(j, PythonRational):
-            j = sympify(j)
+        if mpm:
+            j = sympify(j)._to_mpmath(mp.prec)
         sol_q += j * x0**i
 
     return sol_p / sol_q
 
 
-def _convert_poly_rat(func, x, initcond=True):
+def _convert_poly_rat(func, x, initcond=True, x0=0):
     """Converts Polynomials and Rationals to Holonomic.
     """
 
@@ -1657,6 +1912,10 @@ def _convert_poly_rat(func, x, initcond=True):
     R = QQ.old_poly_ring(x)
     _, Dx = DifferentialOperators(R, 'Dx')
 
+    # if the function is constant
+    if not func.has(x):
+        return HolonomicFunction(Dx, x, 0, func)
+
     if ispoly:
         # differential equation satisfied by polynomial
         sol = func * Dx - func.diff()
@@ -1672,7 +1931,6 @@ def _convert_poly_rat(func, x, initcond=True):
     if not initcond:
         return HolonomicFunction(sol, x)
 
-    x0 = 0
     y0 = _find_conditions(func, x, x0, sol.order)
     while not y0:
         x0 += 1
@@ -1771,6 +2029,8 @@ def _find_conditions(func, x, x0, order):
     y0 = []
     for i in range(order):
         val = func.subs(x, x0)
+        if isinstance(val, NaN):
+            val = limit(func, x, x0)
         if (val.is_finite is not None and not val.is_finite) or isinstance(val, NaN):
             return None
         y0.append(val)
diff --git a/sympy/holonomic/numerical.py b/sympy/holonomic/numerical.py
index c91b6836f5..afbcd1b242 100644
--- a/sympy/holonomic/numerical.py
+++ b/sympy/holonomic/numerical.py
@@ -58,7 +58,7 @@ def _euler(red, x0, x1, y0, a):
     f_0_n = 0
 
     for i in range(a):
-        f_0_n += sympify(DMFsubs(red[i], A))._to_mpmath(mp.prec) * y_0[i]
+        f_0_n += sympify(DMFsubs(red[i], A, mpm=True))._to_mpmath(mp.prec) * y_0[i]
     f_0.append(f_0_n)
 
     sol = []
@@ -85,22 +85,22 @@ def _rk4(red, x0, x1, y0, a):
 
     f_0 = y_0[1:]
     for i in range(a):
-        f_0_n += sympify(DMFsubs(red[i], A))._to_mpmath(mp.prec) * y_0[i]
+        f_0_n += sympify(DMFsubs(red[i], A, mpm=True))._to_mpmath(mp.prec) * y_0[i]
     f_0.append(f_0_n)
 
     f_1 = [y_0[i] + f_0[i]*h/2 for i in range(1, a)]
     for i in range(a):
-        f_1_n += sympify(DMFsubs(red[i], A + h/2))._to_mpmath(mp.prec) * (y_0[i] + f_0[i]*h/2)
+        f_1_n += sympify(DMFsubs(red[i], A + h/2, mpm=True))._to_mpmath(mp.prec) * (y_0[i] + f_0[i]*h/2)
     f_1.append(f_1_n)
 
     f_2 = [y_0[i] + f_1[i]*h/2 for i in range(1, a)]
     for i in range(a):
-        f_2_n += sympify(DMFsubs(red[i], A + h/2))._to_mpmath(mp.prec) * (y_0[i] + f_1[i]*h/2)
+        f_2_n += sympify(DMFsubs(red[i], A + h/2, mpm=True))._to_mpmath(mp.prec) * (y_0[i] + f_1[i]*h/2)
     f_2.append(f_2_n)
 
     f_3 = [y_0[i] + f_2[i]*h for i in range(1, a)]
     for i in range(a):
-        f_3_n += sympify(DMFsubs(red[i], A + h))._to_mpmath(mp.prec) * (y_0[i] + f_2[i]*h)
+        f_3_n += sympify(DMFsubs(red[i], A + h, mpm=True))._to_mpmath(mp.prec) * (y_0[i] + f_2[i]*h)
     f_3.append(f_3_n)
 
     sol = []


Test patch:
diff --git a/sympy/holonomic/tests/test_holonomic.py b/sympy/holonomic/tests/test_holonomic.py
index d917e4c9f7..643df52114 100644
--- a/sympy/holonomic/tests/test_holonomic.py
+++ b/sympy/holonomic/tests/test_holonomic.py
@@ -2,7 +2,7 @@
     DifferentialOperators, from_hyper, from_meijerg, from_sympy)
 from sympy.holonomic.recurrence import RecurrenceOperators, HolonomicSequence
 from sympy import (symbols, hyper, S, sqrt, pi, exp, erf, erfc, sstr,
-    O, I, meijerg, sin, cos, log, cosh, besselj, hyperexpand)
+    O, I, meijerg, sin, cos, log, cosh, besselj, hyperexpand, Ci, EulerGamma, Si)
 from sympy import ZZ, QQ
 
 
@@ -91,8 +91,8 @@ def test_addition_initial_condition():
     assert p + q == r
     p = HolonomicFunction(Dx - x + Dx**2, x, 0, [1, 2])
     q = HolonomicFunction(Dx**2 + x, x, 0, [1, 0])
-    r = HolonomicFunction((-4*x**4 - x**3 - 4*x**2 + 1) + (4*x**3 + x**2 + 3*x + 4)*Dx + \
-        (-6*x + 7)*Dx**2 + (4*x**2 - 7*x + 1)*Dx**3 + (4*x**2 + x + 2)*Dx**4, x, 0, [2, 2, -2, 2])
+    r = HolonomicFunction((-x**4 - x**3/4 - x**2 + 1/4) + (x**3 + x**2/4 + 3*x/4 + 1)*Dx + \
+        (-3*x/2 + 7/4)*Dx**2 + (x**2 - 7*x/4 + 1/4)*Dx**3 + (x**2 + x/4 + 1/2)*Dx**4, x, 0, [2, 2, -2, 2])
     assert p + q == r
     p = HolonomicFunction(Dx**2 + 4*x*Dx + x**2, x, 0, [3, 4])
     q = HolonomicFunction(Dx**2 + 1, x, 0, [1, 1])
@@ -105,6 +105,11 @@ def test_addition_initial_condition():
     r = HolonomicFunction((-x**2 - x + 1) + (x**2 + x)*Dx + (-x - 2)*Dx**3 + \
         (x + 1)*Dx**4, x, 2, [4, 1, 2, -5 ])
     assert p + q == r
+    p = from_sympy(sin(x))
+    q = from_sympy(1/x)
+    r = HolonomicFunction((x**2 + 6) + (x**3 + 2*x)*Dx + (x**2 + 6)*Dx**2 + (x**3 + 2*x)*Dx**3, \
+        x, 1, [sin(1) + 1, -1 + cos(1), -sin(1) + 2])
+    assert p + q == r
 
 def test_multiplication_initial_condition():
     x = symbols('x')
@@ -116,12 +121,12 @@ def test_multiplication_initial_condition():
     assert p * q == r
     p = HolonomicFunction(Dx**2 + x, x, 0, [1, 0])
     q = HolonomicFunction(Dx**3 - x**2, x, 0, [3, 3, 3])
-    r = HolonomicFunction((27*x**8 - 37*x**7 - 10*x**6 - 492*x**5 - 552*x**4 + 160*x**3 + \
-        1212*x**2 + 216*x + 360) + (162*x**7 - 384*x**6 - 294*x**5 - 84*x**4 + 24*x**3 + \
-        756*x**2 + 120*x - 1080)*Dx + (81*x**6 - 246*x**5 + 228*x**4 + 36*x**3 + \
-        660*x**2 - 720*x)*Dx**2 + (-54*x**6 + 128*x**5 - 18*x**4 - 240*x**2 + 600)*Dx**3 + \
-        (81*x**5 - 192*x**4 - 84*x**3 + 162*x**2 - 60*x - 180)*Dx**4 + (-108*x**3 + \
-        192*x**2 + 72*x)*Dx**5 + (27*x**4 - 64*x**3 - 36*x**2 + 60)*Dx**6, x, 0, [3, 3, 3, -3, -12, -24])
+    r = HolonomicFunction((x**8 - 37*x**7/27 - 10*x**6/27 - 164*x**5/9 - 184*x**4/9 + \
+        160*x**3/27 + 404*x**2/9 + 8*x + 40/3) + (6*x**7 - 128*x**6/9 - 98*x**5/9 - 28*x**4/9 + \
+        8*x**3/9 + 28*x**2 + 40*x/9 - 40)*Dx + (3*x**6 - 82*x**5/9 + 76*x**4/9 + 4*x**3/3 + \
+        220*x**2/9 - 80*x/3)*Dx**2 + (-2*x**6 + 128*x**5/27 - 2*x**4/3 -80*x**2/9 + 200/9)*Dx**3 + \
+        (3*x**5 - 64*x**4/9 - 28*x**3/9 + 6*x**2 - 20*x/9 - 20/3)*Dx**4 + (-4*x**3 + 64*x**2/9 + \
+            8*x/3)*Dx**5 + (x**4 - 64*x**3/27 - 4*x**2/3 + 20/9)*Dx**6, x, 0, [3, 3, 3, -3, -12, -24])
     assert p * q == r
     p = HolonomicFunction(Dx - 1, x, 0, [2])
     q = HolonomicFunction(Dx**2 + 1, x, 0, [0, 1])
@@ -134,6 +139,10 @@ def test_multiplication_initial_condition():
     q = HolonomicFunction(Dx**3 + 1, x, 0, [1, 2, 1])
     r = HolonomicFunction(6*Dx + 3*Dx**2 + 2*Dx**3 - 3*Dx**4 + Dx**6, x, 0, [1, 5, 14, 17, 17, 2])
     assert p * q == r
+    p = from_sympy(sin(x))
+    q = from_sympy(1/x)
+    r = HolonomicFunction(x + 2*Dx + x*Dx**2, x, 1, [sin(1), -sin(1) + cos(1)])
+    assert p * q == r
 
 def test_HolonomicFunction_composition():
     x = symbols('x')
@@ -168,7 +177,7 @@ def test_from_hyper():
     r = from_hyper(p)
     assert r == q
     p = from_hyper(hyper([1], [S(3)/2], x**2/4))
-    q = HolonomicFunction(-2*x + (-x**2 + 4)*Dx + 2*x*Dx**2, x)
+    q = HolonomicFunction(-x + (-x**2/2 + 2)*Dx + x*Dx**2, x)
     x0 = 1
     y0 = '[sqrt(pi)*exp(1/4)*erf(1/2), -sqrt(pi)*exp(1/4)*erf(1/2)/2 + 1]'
     assert sstr(p.y0) == y0
@@ -178,14 +187,14 @@ def test_from_meijerg():
     x = symbols('x')
     R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
     p = from_meijerg(meijerg(([], [S(3)/2]), ([S(1)/2], [S(1)/2, 1]), x))
-    q = HolonomicFunction(2*x - 1 + (-4*x**2 + x)*Dx + 4*x**2*Dx**2 + 4*x**3*Dx**3, x, 1, \
+    q = HolonomicFunction(x/2 - 1/4 + (-x**2 + x/4)*Dx + x**2*Dx**2 + x**3*Dx**3, x, 1, \
         [1/sqrt(pi), 1/(2*sqrt(pi)), -1/(4*sqrt(pi))])
     assert p == q
     p = from_meijerg(meijerg(([], []), ([0], []), x))
     q = HolonomicFunction(1 + Dx, x, 0, 1)
     assert p == q
     p = from_meijerg(meijerg(([1], []), ([S(1)/2], [0]), x))
-    q = HolonomicFunction((2*x + 1)*Dx + 2*x*Dx**2, x, 1, [sqrt(pi)*erf(1), exp(-1)])
+    q = HolonomicFunction((x + 1/2)*Dx + x*Dx**2, x, 1, [sqrt(pi)*erf(1), exp(-1)])
     assert p == q
     p = from_meijerg(meijerg(([0], [1]), ([0], []), 2*x**2))
     q = HolonomicFunction((3*x**2 - 1)*Dx + x**3*Dx**2, x, 1, [-exp(-S(1)/2) + 1, -exp(-S(1)/2)])
@@ -200,13 +209,13 @@ def test_to_Sequence():
     q = (HolonomicSequence(1 + (n + 2)*Sn**2 + (n**4 + 6*n**3 + 11*n**2 + 6*n)*Sn**3), 1)
     assert p == q
     p = HolonomicFunction(x**2*Dx**4 + x**3 + Dx**2, x).to_sequence()
-    q = (HolonomicSequence(1 + (n**4 + 14*n**3 + 72*n**2 + 163*n + 140)*Sn**5, n), 0)
+    q = (HolonomicSequence(1 + (n**4 + 14*n**3 + 72*n**2 + 163*n + 140)*Sn**5), 0)
     assert p == q
     p = HolonomicFunction(x**3*Dx**4 + 1 + Dx**2, x).to_sequence()
-    q = (HolonomicSequence(1 + (n**4 - 2*n**3 - n**2 + 2*n)*Sn + (n**2 + 3*n + 2)*Sn**2, n), 3)
+    q = (HolonomicSequence(1 + (n**4 - 2*n**3 - n**2 + 2*n)*Sn + (n**2 + 3*n + 2)*Sn**2), 3)
     assert p == q
     p = HolonomicFunction(3*x**3*Dx**4 + 2*x*Dx + x*Dx**3, x).to_sequence()
-    q = (HolonomicSequence(2*n + (3*n**4 - 6*n**3 - 3*n**2 + 6*n)*Sn + (n**3 + 3*n**2 + 2*n)*Sn**2, n), 3)
+    q = (HolonomicSequence(2*n + (3*n**4 - 6*n**3 - 3*n**2 + 6*n)*Sn + (n**3 + 3*n**2 + 2*n)*Sn**2), 3)
     assert p == q
 
 def test_to_Sequence_Initial_Coniditons():
@@ -226,6 +235,16 @@ def test_to_Sequence_Initial_Coniditons():
     p = HolonomicFunction(x**3*Dx**5 + 1 + Dx, x).to_sequence()
     q = (HolonomicSequence(1 + (n + 1)*Sn + (n**5 - 5*n**3 + 4*n)*Sn**2), 3)
     assert p == q
+    C_1, C_2, C_3 = symbols('C_1, C_2, C_3')
+    p = from_sympy(log(1+x**2))
+    q = (HolonomicSequence(n**2 + (n**2 + 2*n)*Sn**2, [0, 0, C_2, 0]), 2)
+    assert p.to_sequence() == q
+    p = p.diff()
+    q = (HolonomicSequence((n + 1) + (n + 1)*Sn**2, [0, C_1, 0]), 1)
+    assert p.to_sequence() == q
+    p = from_sympy(erf(x) + x).to_sequence()
+    q = (HolonomicSequence((2*n**2 - 2*n) + (n**3 + 2*n**2 - n - 2)*Sn**2, [0, 1 + 2/sqrt(pi), 0, C_3]), 2)
+    assert p == q
 
 def test_series():
     x = symbols('x')
@@ -251,6 +270,10 @@ def test_series():
         (4-6*x**3+2*x**4)*Dx**2, x, 0, [1, 0]).series(n=7)
     q = 1 - 3*x**2/4 - x**3/4 - 5*x**4/32 - 3*x**5/40 - 17*x**6/384 + O(x**7)
     assert p == q
+    p = from_sympy(erf(x) + x).series(n=10)
+    C_3 = symbols('C_3')
+    q = (erf(x) + x).series(n=10)
+    assert p.subs(C_3, -2/(3*sqrt(pi))) == q
 
 def test_evalf_euler():
     x = symbols('x')
@@ -410,17 +433,17 @@ def test_from_sympy():
     x = symbols('x')
     R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
     p = from_sympy((sin(x)/x)**2)
-    q = HolonomicFunction(8*x + (4*x**2 + 6)*Dx + 6*x*Dx**2 + x**2*Dx**3, x, 1, \
-        [sin(1)**2, -2*sin(1)**2 + 2*sin(1)*cos(1), -8*sin(1)*cos(1) + 2*cos(1)**2 + 4*sin(1)**2])
+    q = HolonomicFunction(8*x + (4*x**2 + 6)*Dx + 6*x*Dx**2 + x**2*Dx**3, x, 0, \
+        [1, 0, -2/3])
     assert p == q
     p = from_sympy(1/(1+x**2)**2)
     q = HolonomicFunction(4*x + (x**2 + 1)*Dx, x, 0, 1)
     assert p == q
     p = from_sympy(exp(x)*sin(x)+x*log(1+x))
-    q = HolonomicFunction((4*x**3 + 20*x**2 + 40*x + 36) + (-4*x**4 - 20*x**3 - 40*x**2 \
-        - 36*x)*Dx + (4*x**5 + 12*x**4 + 14*x**3 + 16*x**2 + 20*x - 8)*Dx**2 + \
-        (-4*x**5 - 10*x**4 - 4*x**3 + 4*x**2 - 2*x + 8)*Dx**3 + (2*x**5 + 4*x**4 - 2*x**3 - \
-        7*x**2 + 2*x + 5)*Dx**4, x, 0, [0, 1, 4, -1])
+    q = HolonomicFunction((2*x**3 + 10*x**2 + 20*x + 18) + (-2*x**4 - 10*x**3 - 20*x**2 \
+        - 18*x)*Dx + (2*x**5 + 6*x**4 + 7*x**3 + 8*x**2 + 10*x - 4)*Dx**2 + \
+        (-2*x**5 - 5*x**4 - 2*x**3 + 2*x**2 - x + 4)*Dx**3 + (x**5 + 2*x**4 - x**3 - \
+        7*x**2/2 + x + 5/2)*Dx**4, x, 0, [0, 1, 4, -1])
     assert p == q
     p = from_sympy(x*exp(x)+cos(x)+1)
     q = HolonomicFunction((-x - 3)*Dx + (x + 2)*Dx**2 + (-x - 3)*Dx**3 + (x + 2)*Dx**4, x, \
@@ -431,8 +454,8 @@ def test_from_sympy():
     q = HolonomicFunction(Dx + (3*x + 3)*Dx**2 + (x**2 + 2*x + 1)*Dx**3, x, 0, [1, 0, 2])
     assert p == q
     p = from_sympy(erf(x)**2 + x)
-    q = HolonomicFunction((32*x**4 - 8*x**2 + 8)*Dx**2 + (24*x**3 - 2*x)*Dx**3 + \
-        (4*x**2+ 1)*Dx**4, x, 0, [0, 1, 8/pi, 0])
+    q = HolonomicFunction((8*x**4 - 2*x**2 + 2)*Dx**2 + (6*x**3 - x/2)*Dx**3 + \
+        (x**2+ 1/4)*Dx**4, x, 0, [0, 1, 8/pi, 0])
     assert p == q
     p = from_sympy(cosh(x)*x)
     q = HolonomicFunction((-x**2 + 2) -2*x*Dx + x**2*Dx**2, x, 0, [0, 1])
@@ -441,9 +464,24 @@ def test_from_sympy():
     q = HolonomicFunction((x**2 - 4) + x*Dx + x**2*Dx**2, x, 0, [0, 0])
     assert p == q
     p = from_sympy(besselj(0, x) + exp(x))
-    q = HolonomicFunction((-2*x**2 - x + 1) + (2*x**2 - x - 3)*Dx + (-2*x**2 + x + 2)*Dx**2 +\
-        (2*x**2 + x)*Dx**3, x, 0, [2, 1, 1/2])
+    q = HolonomicFunction((-x**2 - x/2 + 1/2) + (x**2 - x/2 - 3/2)*Dx + (-x**2 + x/2 + 1)*Dx**2 +\
+        (x**2 + x/2)*Dx**3, x, 0, [2, 1, 1/2])
+    assert p == q
+    p = from_sympy(sin(x)**2/x)
+    q = HolonomicFunction(4 + 4*x*Dx + 3*Dx**2 + x*Dx**3, x, 0, [0, 1, 0])
     assert p == q
+    p = from_sympy(sin(x)**2/x, x0=2)
+    q = HolonomicFunction((4) + (4*x)*Dx + (3)*Dx**2 + (x)*Dx**3, x, 2, [sin(2)**2/2,
+        sin(2)*cos(2) - sin(2)**2/4, -3*sin(2)**2/4 + cos(2)**2 - sin(2)*cos(2)])
+    assert p == q
+    p = from_sympy(log(x)/2 - Ci(2*x)/2 + Ci(2)/2)
+    q = HolonomicFunction(4*Dx + 4*x*Dx**2 + 3*Dx**3 + x*Dx**4, x, 0, \
+        [-log(2)/2 - EulerGamma/2 + Ci(2)/2, 0, 1, 0])
+    assert p == q
+    p = p.to_sympy()
+    q = log(x)/2 - Ci(2*x)/2 + Ci(2)/2
+    assert p == q
+
 
 def test_to_hyper():
     x = symbols('x')
@@ -487,3 +525,51 @@ def test_to_sympy():
         (x**2 - x)*Dx**2, x, 0, [1, 2]).to_sympy()
     q = 1/(x**2 - 2*x + 1)
     assert p == q
+    p = from_sympy(sin(x)**2/x).integrate((x, 0, x)).to_sympy()
+    q = (sin(x)**2/x).integrate((x, 0, x))
+    assert p == q
+    C_1, C_2, C_3 = symbols('C_1, C_2, C_3')
+    p = from_sympy(log(1+x**2)).to_sympy()
+    q = C_2*log(x**2 + 1)
+    assert p == q
+    p = from_sympy(log(1+x**2)).diff().to_sympy()
+    q = C_1*x/(x**2 + 1)
+    assert p == q
+    p = from_sympy(erf(x) + x).to_sympy()
+    q = 3*C_3*x - 3*sqrt(pi)*C_3*erf(x)/2 + x + 2*x/sqrt(pi)
+    assert p == q
+
+def test_integrate():
+    x = symbols('x')
+    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
+    p = from_sympy(sin(x)**2/x, x0=1).integrate((x, 2, 3))
+    q = '0.166270406994788'
+    assert sstr(p) == q
+    p = from_sympy(sin(x)).integrate((x, 0, x)).to_sympy()
+    q = 1 - cos(x)
+    assert p == q
+    p = from_sympy(sin(x)).integrate((x, 0, 3))
+    q = '1.98999246812687'
+    assert sstr(p) == q
+    p = from_sympy(sin(x)/x, x0=1).integrate((x, 1, 2))
+    q = '0.659329913368450'
+    assert sstr(p) == q
+    p = from_sympy(sin(x)**2/x, x0=1).integrate((x, 1, 0))
+    q = '-0.423690480850035'
+    assert sstr(p) == q
+
+def test_diff():
+    x, y = symbols('x, y')
+    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
+    p = HolonomicFunction(x*Dx**2 + 1, x, 0, [0, 1])
+    assert p.diff().to_sympy() == p.to_sympy().diff().simplify()
+    p = HolonomicFunction(Dx**2 - 1, x, 0, [1, 0])
+    assert p.diff(x, 2).to_sympy() == p.to_sympy()
+    p = from_sympy(Si(x))
+    assert p.diff().to_sympy() == sin(x)/x
+    assert p.diff(y) == 0
+    C_1, C_2, C_3 = symbols('C_1, C_2, C_3')
+    q = Si(x)
+    assert p.diff(x).to_sympy() == q.diff()
+    assert p.diff(x, 2).to_sympy().subs(C_1, -S(1)/3) == q.diff(x, 2).simplify()
+    assert p.diff(x, 3).series().subs(C_2, S(1)/10) == q.diff(x, 3).series()