What does r() function mean in the return value of SymPy's dsolve?

Question

I want to evaluate the value of phi(+oo) where phi(xi) is the solution of ODE

Eq(Derivative(phi(xi), (xi, 2)), (-K + xi**2)*phi(xi))

and K is a known real variable. By dsolve, I got the solution:

Eq(phi(xi), -K*xi**5*r(3)/20 + C2*(K**2*xi**4/24 - K*xi**2/2 + xi**4/12 + 1) + C1*xi*(xi**4/20 + 1) + O(xi**6))

with an unknown function r() in the first term on the right-hand side. Here is my code:

import numpy as np
import matplotlib.pyplot as plt
import sympy
from sympy import I, pi, oo

sympy.init_printing()

def apply_ics(sol, ics, x, known_params):
    """
    Apply the initial conditions (ics), given as a dictionary on
    the form ics = {y(0): y0, y(x).diff(x).subs(x, 0): yp0, ...},
    to the solution of the ODE with independent variable x.
    The undetermined integration constants C1, C2, ... are extracted
    from the free symbols of the ODE solution, excluding symbols in
    the known_params list.
    """

    free_params = sol.free_symbols - set(known_params)
    eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics)
            for n in range(len(ics))]
    sol_params = sympy.solve(eqs, free_params)

    return sol.subs(sol_params)

K = sympy.Symbol('K', positive = True)
xi = sympy.Symbol('xi',real = True)
phi = sympy.Function('phi')
ode = sympy.Eq( phi(xi).diff(xi, 2), (xi**2-K)*phi(xi))

ode_sol = sympy.dsolve(ode)
ics = { phi(0):1, phi(xi).diff(xi).subs(xi,0): 0}
phi_xi_sol = apply_ics(ode_sol, ics, xi, [K])

Where ode_sol is the solution, phi_xi_sol is the solution after initial conditions are applied. Since r() is undefined in NumPy I can't evaluate the results by

for g in [0.9, 0.95, 1, 1.05, 1.2]:
    phi_xi = sympy.lambdify(xi, phi_xi_sol.rhs.subs({K:g}), 'numpy')

Does anyone know what this function r() mean and how should I deal with it?

Answer 1

As visible in the form of the result, the solver falls back to a power series solution (instead of searching the solution in terms of parabolic cylinder functions as WolframAlpha does).

So let's set phi(xi)=sum a[k]*xi^k leading to the coefficient equations (using a[k]=0 for k<0)

(k+2)(k+1)a[k+2] = -K*a[k] + a[k-2]

a[0] = C2  
a[1] = C1  
a[2] = -K/2*C2
a[3] = -K/6*C1
a[4] = (K^2/2 + 1)/12*C2
a[5] = (K^2/6 + 1)/20*C1

Inserting that the the power series solution should have been

C2*(1-K/2*xi**2+(K**2/24+1/12)*xi**4) + C1*xi*(1-K/6*xi**2+(K/120+1/20)*xi**4) + O(xi**6)

Comparing with the sympy solution, all terms containing both C1 and K are missing, especially the missing degree 3 term is not explainable. It seems that the solution process was prematurely ended, or some equation transformation was not correctly reversed.

Please note that the ODE solver routines in sympy are experimental and rudimentary. Also, the power series solution gives only valid information for small values of xi, there is no way to derive any exact value for the limit at +oo.