code: def solve_aux_eq(numa, dena, numy, deny, x, m):
      """
      Helper function to find a polynomial solution
      of degree m for the auxiliary differential
      equation.
      """
      # Assume that the solution is of the type
      # p(x) = C_0 + C_1*x + ... + C_{m-1}*x**(m-1) + x**m
      psyms = symbols(f"C0:{m}", cls=Dummy)
      K = ZZ[psyms]
      psol = Poly(K.gens, x, domain=K) + Poly(x**m, x, domain=K)
  
      # Eq (5.16) in Thesis - Pg 81
      auxeq = (
          dena * (numy.diff(x) * deny - numy * deny.diff(x) + numy**2)
          - numa * deny**2
      ) * psol
      if m >= 1:
          px = psol.diff(x)
          auxeq += px * (2 * numy * deny * dena)
      if m >= 2:
          auxeq += px.diff(x) * (deny**2 * dena)
      if m != 0:
          # m is a non-zero integer. Find the constant terms using undetermined coefficients
          return psol, linsolve_dict(auxeq.all_coeffs(), psyms), True
      else:
          # m == 0 . Check if 1 (x**0) is a solution to the auxiliary equation
          return S.One, auxeq, auxeq == 0
  
imports: from itertools import product
  from sympy.core import S
  from sympy.core.add import Add
  from sympy.core.numbers import oo, Float
  from sympy.core.function import count_ops
  from sympy.core.relational import Eq
  from sympy.core.symbol import symbols, Symbol, Dummy
  from sympy.functions import sqrt, exp
  from sympy.functions.elementary.complexes import sign
  from sympy.integrals.integrals import Integral
  from sympy.polys.domains import ZZ
  from sympy.polys.polytools import Poly
  from sympy.polys.polyroots import roots
  from sympy.solvers.solveset import linsolve
label_desc: Use image node_id 20 for calling a global function with example usage: solve_aux_eq(numa, dena, numy, deny, x, m) and returns: psol,  linsolve_dict,  True, S,  auxeq,  unknown