'Returns function after composition of a holonomic
function with an algebraic function. The method can\'t compute
initial conditions for the result by itself, so they can be also be
provided.
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).composition(x**2, 0, [1])  # e^(x**2)
HolonomicFunction((-2*x) + (1)*Dx, x, 0, [1])
>>> HolonomicFunction(Dx**2 + 1, x).composition(x**2 - 1, 1, [1, 0])
HolonomicFunction((4*x**3) + (-1)*Dx + (x)*Dx**2, x, 1, [1, 0])
See Also
from_hyper()'

def composition(self, expr, *args, **kwargs):

    R = self.annihilator.parent
    a = self.annihilator.order
    diff = expr.diff(self.x)
    listofpoly = self.annihilator.listofpoly
    for (i, j) in enumerate(listofpoly):
        if isinstance(j, self.annihilator.parent.base.dtype):
            listofpoly[i] = self.annihilator.parent.base.to_sympy(j)
    r = listofpoly[a].subs({self.x: expr})
    subs = [((- listofpoly[i].subs({self.x: expr})) / r) for i in range(a)]
    coeffs = [S(0) for i in range(a)]
    coeffs[0] = S(1)
    system = [coeffs]
    homogeneous = Matrix([[S(0) for i in range(a)]]).transpose()
    sol = S(0)
    while sol.is_zero:
        coeffs_next = [p.diff(self.x) for p in coeffs]
        for i in range((a - 1)):
            coeffs_next[(i + 1)] += (coeffs[i] * diff)
        for i in range(a):
            coeffs_next[i] += ((coeffs[(-1)] * subs[i]) * diff)
        coeffs = coeffs_next
        system.append(coeffs)
        sol_tuple = Matrix(system).transpose().gauss_jordan_solve(homogeneous)
        sol = sol_tuple[0]
    tau = sol.atoms(Dummy).pop()
    sol = sol.subs(tau, 1)
    sol = _normalize(sol[0:], R, negative=False)
    if args:
        return HolonomicFunction(sol, self.x, args[0], args[1])
    return HolonomicFunction(sol, self.x)
