'Finds the power series expansion of given holonomic function about :math:`x_0`.
A list of series might be returned if :math:`x_0` is a regular point with
multiple roots of the indicial equation.
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\')
>>> HolonomicFunction(Dx - 1, x, 0, [1]).series()  # e^x
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
>>> HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]).series(n=8)  # sin(x)
x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)
See Also
HolonomicFunction.to_sequence()'

def series(self, n=6, coefficient=False, order=True, _recur=None):

    if (_recur == None):
        recurrence = self.to_sequence()
    else:
        recurrence = _recur
    if (isinstance(recurrence, tuple) and (len(recurrence) == 2)):
        recurrence = recurrence[0]
        constantpower = 0
    elif (isinstance(recurrence, tuple) and (len(recurrence) == 3)):
        constantpower = recurrence[1]
        recurrence = recurrence[0]
    elif ((len(recurrence) == 1) and (len(recurrence[0]) == 2)):
        recurrence = recurrence[0][0]
        constantpower = 0
    elif ((len(recurrence) == 1) and (len(recurrence[0]) == 3)):
        constantpower = recurrence[0][1]
        recurrence = recurrence[0][0]
    else:
        sol = []
        for i in recurrence:
            sol.append(self.series(_recur=i))
        return sol
    n = (n - int(constantpower))
    l = (len(recurrence.u0) - 1)
    k = recurrence.recurrence.order
    x = self.x
    x0 = self.x0
    seq_dmp = recurrence.recurrence.listofpoly
    R = recurrence.recurrence.parent.base
    K = R.get_field()
    seq = []
    for (i, j) in enumerate(seq_dmp):
        seq.append(K.new(j.rep))
    sub = [((- seq[i]) / seq[k]) for i in range(k)]
    sol = [i for i in recurrence.u0]
    if ((l + 1) >= n):
        pass
    else:
        for i in range(((l + 1) - k), (n - k)):
            coeff = S(0)
            for j in range(k):
                if ((i + j) >= 0):
                    coeff += (DMFsubs(sub[j], i) * sol[(i + j)])
            sol.append(coeff)
    if coefficient:
        return sol
    ser = S(0)
    for (i, j) in enumerate(sol):
        ser += ((x ** (i + constantpower)) * j)
    if order:
        ser += Order((x ** (n + int(constantpower))), x)
    if (x0 != 0):
        return ser.subs(x, (x - x0))
    return ser
