code: def rsolve(f, y, init=None):
      r"""
      Solve univariate recurrence with rational coefficients.
  
      Given `k`-th order linear recurrence `\operatorname{L} y = f`,
      or equivalently:
  
      .. math:: a_{k}(n) y(n+k) + a_{k-1}(n) y(n+k-1) +
                \cdots + a_{0}(n) y(n) = f(n)
  
      where `a_{i}(n)`, for `i=0, \ldots, k`, are polynomials or rational
      functions in `n`, and `f` is a hypergeometric function or a sum
      of a fixed number of pairwise dissimilar hypergeometric terms in
      `n`, finds all solutions or returns ``None``, if none were found.
  
      Initial conditions can be given as a dictionary in two forms:
  
          (1) ``{  n_0  : v_0,   n_1  : v_1, ...,   n_m  : v_m}``
          (2) ``{y(n_0) : v_0, y(n_1) : v_1, ..., y(n_m) : v_m}``
  
      or as a list ``L`` of values:
  
          ``L = [v_0, v_1, ..., v_m]``
  
      where ``L[i] = v_i``, for `i=0, \ldots, m`, maps to `y(n_i)`.
  
      Examples
      ========
  
      Lets consider the following recurrence:
  
      .. math:: (n - 1) y(n + 2) - (n^2 + 3 n - 2) y(n + 1) +
                2 n (n + 1) y(n) = 0
  
      >>> from sympy import Function, rsolve
      >>> from sympy.abc import n
      >>> y = Function('y')
  
      >>> f = (n - 1)*y(n + 2) - (n**2 + 3*n - 2)*y(n + 1) + 2*n*(n + 1)*y(n)
  
      >>> rsolve(f, y(n))
      2**n*C0 + C1*factorial(n)
  
      >>> rsolve(f, y(n), {y(0):0, y(1):3})
      3*2**n - 3*factorial(n)
  
      See Also
      ========
  
      rsolve_poly, rsolve_ratio, rsolve_hyper
  
      """
      if isinstance(f, Equality):
          f = f.lhs - f.rhs
  
      n = y.args[0]
      k = Wild("k", exclude=(n,))
  
      # Preprocess user input to allow things like
      # y(n) + a*(y(n + 1) + y(n - 1))/2
      f = f.expand().collect(y.func(Wild("m", integer=True)))
  
      h_part = defaultdict(list)
      i_part = []
      for g in Add.make_args(f):
          coeff, dep = g.as_coeff_mul(y.func)
          if not dep:
              i_part.append(coeff)
              continue
          for h in dep:
              if h.is_Function and h.func == y.func:
                  result = h.args[0].match(n + k)
                  if result is not None:
                      h_part[int(result[k])].append(coeff)
                      continue
              raise ValueError(
                  "'%s(%s + k)' expected, got '%s'" % (y.func, n, h)
              )
      for k in h_part:
          h_part[k] = Add(*h_part[k])
      h_part.default_factory = lambda: 0
      i_part = Add(*i_part)
  
      for k, coeff in h_part.items():
          h_part[k] = simplify(coeff)
  
      common = S.One
  
      if (
          not i_part.is_zero
          and not i_part.is_hypergeometric(n)
          and not (
              i_part.is_Add
              and all(
                  (x.is_hypergeometric(n) for x in i_part.expand().args)
              )
          )
      ):
          raise ValueError(
              "The independent term should be a sum of hypergeometric functions, got '%s'"
              % i_part
          )
  
      for coeff in h_part.values():
          if coeff.is_rational_function(n):
              if not coeff.is_polynomial(n):
                  common = lcm(common, coeff.as_numer_denom()[1], n)
          else:
              raise ValueError(
                  "Polynomial or rational function expected, got '%s'"
                  % coeff
              )
  
      i_numer, i_denom = i_part.as_numer_denom()
  
      if i_denom.is_polynomial(n):
          common = lcm(common, i_denom, n)
  
      if common is not S.One:
          for k, coeff in h_part.items():
              numer, denom = coeff.as_numer_denom()
              h_part[k] = numer * quo(common, denom, n)
  
          i_part = i_numer * quo(common, i_denom, n)
  
      K_min = min(h_part.keys())
  
      if K_min < 0:
          K = abs(K_min)
  
          H_part = defaultdict(lambda: S.Zero)
          i_part = i_part.subs(n, n + K).expand()
          common = common.subs(n, n + K).expand()
  
          for k, coeff in h_part.items():
              H_part[k + K] = coeff.subs(n, n + K).expand()
      else:
          H_part = h_part
  
      K_max = max(H_part.keys())
      coeffs = [H_part[i] for i in range(K_max + 1)]
  
      result = rsolve_hyper(coeffs, -i_part, n, symbols=True)
  
      if result is None:
          return None
  
      solution, symbols = result
  
      if init in ({}, []):
          init = None
  
      if symbols and init is not None:
          if isinstance(init, list):
              init = {i: init[i] for i in range(len(init))}
  
          equations = []
  
          for k, v in init.items():
              try:
                  i = int(k)
              except TypeError:
                  if k.is_Function and k.func == y.func:
                      i = int(k.args[0])
                  else:
                      raise ValueError(
                          "Integer or term expected, got '%s'" % k
                      )
  
              eq = solution.subs(n, i) - v
              if eq.has(S.NaN):
                  eq = solution.limit(n, i) - v
              equations.append(eq)
  
          result = solve(equations, *symbols)
  
          if not result:
              return None
          else:
              solution = solution.subs(result)
  
      return solution
  
imports: from collections import defaultdict
  from sympy.concrete import product
  from sympy.core.singleton import S
  from sympy.core.numbers import Rational, I
  from sympy.core.symbol import Symbol, Wild, Dummy
  from sympy.core.relational import Equality
  from sympy.core.add import Add
  from sympy.core.mul import Mul
  from sympy.core.sorting import default_sort_key
  from sympy.core.sympify import sympify
  from sympy.simplify import simplify, hypersimp, hypersimilar
  from sympy.solvers import solve, solve_undetermined_coeffs
  from sympy.polys import Poly, quo, gcd, lcm, roots, resultant
  from sympy.functions import binomial, factorial, FallingFactorial, RisingFactorial
  from sympy.matrices import Matrix, casoratian
  from sympy.utilities.iterables import numbered_symbols
label_desc: Use image node_id 4 for calling a global function with example usage: rsolve(f, y, init) and returns: solution, None, None