Repository: sympy/sympy
Pull request: Add solvers for Bessel and Airy differential equations
#### References to other Issues or PRs

closes #10870
fixes #9819

#### Brief description of what is fixed or changed
Currently, sympy does not support the ODEs solution in term of Bessel function.
So for this, I have taken help from pr #10870 and Implemented  ODEs in term of Bessel and Airy's function.

```
Eq(x ** 2*Derivative(y(x),x,x)-(ax+b)*y(x))
The above equation now has solution c1*Ai(x)+x2*Bi(x)
where Ai(x) and Bi(x) are Airy's functions of first and second kind respectively

```
```
Eq(x**2*f(x).diff(x,2)+x*f(x).diff(x)+(x**2-a**2)*f(x),0)
The above equation now has solution 
                 ⎛   ____   ⎞             ⎛   ____   ⎞
                 ⎜  ╱  2    ⎟             ⎜  ╱  2    ⎟
f(x) = C₁⋅besselj⎝╲╱  a  , x⎠ + C₂⋅bessely⎝╲╱  a  , x⎠


where besselj(n,x) and bessely(n,x) are Bessel's functions of first and second kind respectively
```

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- solvers
  - Add solvers for Bessel and Airy differential equations in dsolve

<!-- END RELEASE NOTES -->


Problem statement:
ODE solver cannot solve Bessel equation
I tried to solve the Bessel equation with Sympy (version 0.7.6), but got a weird result (the Bessel function is expected):

```
>>> x = symbols('x')
>>> f = Function('f')
>>> eq=Derivative(f(x),x,2)*x*x+x*Derivative(f(x),x,1)+(x*x)*f(x)
>>> dsolve(eq)
f(x) == C1*(x**4/64 - x**2/4 + 1) + O(x**6)
```

Hints from issue discussion:
Yes, we need somebody to implement a routine in ode.py that will recognize and solve the 
classical 2nd-order variable-coefficient linear ordinary differential equations like the Bessel equation and the Airy equation.  The solution sympy is currently giving is the series for Bessel's J function of order 0:

```
In [51]: besselj(0,x).series(x,0,6)
Out[51]: 
     2    4        
    x    x     ⎛ 6⎞
1 - ── + ── + O⎝x ⎠
    4    64       
```

The full solution should be `C1*bessely(0,x)+C2*besselj(0,x)`.

I would like to work on this.

@nagarakshay  I [started to work](https://github.com/zanzibar7/sympy.git) on this several months ago and got a little way, but not to Bessel functions.  Airy equations seemed like an easier place to start, and then Bessel equations could be added on.  Ultimately, hypergeometric equations could also be addressed.  Please feel free to look at my changes, and feel free to incorporate them.

@zanzibar7 Thank you. I will surely have a look and incorporate them in my PR.

@zanzibar7 What further changes may be done?

See PR #10870

Commit message: add checkodesol in test case and some code cleaning

Patch to review:
diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index d00b13ad520f..f526b93ff3cb 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -235,7 +235,7 @@
 from collections import defaultdict
 from itertools import islice
 
-from sympy.core import Add, S, Mul, Pow, oo
+from sympy.core import Add, S, Mul, Pow, oo, Rational
 from sympy.core.compatibility import ordered, iterable, is_sequence, range, string_types
 from sympy.core.containers import Tuple
 from sympy.core.exprtools import factor_terms
@@ -251,7 +251,7 @@
 from sympy.logic.boolalg import (BooleanAtom, And, Not, BooleanTrue,
                                 BooleanFalse)
 from sympy.functions import cos, exp, im, log, re, sin, tan, sqrt, \
-    atan2, conjugate, Piecewise
+    atan2, conjugate, Piecewise, cbrt, besselj, bessely, airyai, airybi
 from sympy.functions.combinatorial.factorials import factorial
 from sympy.integrals.integrals import Integral, integrate
 from sympy.matrices import wronskian, Matrix, eye, zeros
@@ -307,6 +307,9 @@
     "nth_linear_constant_coeff_variation_of_parameters",
     "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters",
     "Liouville",
+    "order_reducible",
+    "2nd_linear_airy",
+    "2nd_linear_bessel",
     "nth_order_reducible",
     "2nd_power_series_ordinary",
     "2nd_power_series_regular",
@@ -638,7 +641,6 @@ def recur_len(l):
         hints = _desolve(eq, func=func,
             hint=hint, simplify=True, xi=xi, eta=eta, type='ode', ics=ics,
             x0=x0, n=n, **kwargs)
-
         eq = hints.pop('eq', eq)
         all_ = hints.pop('all', False)
         if all_:
@@ -1351,8 +1353,39 @@ class in it.  Note that a hint may do this anyway if
                             not check.has(zoo) and not check.has(-oo):
                             coeff_dict = {'p': p, 'q': q, 'x0': point, 'terms': terms}
                             matching_hints["2nd_power_series_regular"] = coeff_dict
-
-
+                            # If the ODE has regular singular point at x0 and is of the form
+                            # Eq((x)**2*Derivative(y(x),x,x)+x*Derivative(y(x),x)+(a4**2*x**2-n**2)*y(x) thus Bessel's equation
+                            if (p.is_constant(x) and p != 0) and r[c3] != 0:
+                                a4 = Wild('a4', exclude=[x,f(x),df])
+                                b4 = Wild('b4', exclude=[x,f(x),df])
+                                c4 = Wild('b4', exclude=[x,f(x),df])
+                                coeff = r[a3].match(a4*x**2)
+                                if coeff:
+                                    coeff = coeff[a4]
+                                    if not coeff.has(x):
+                                        r[c3] = r[c3]/coeff
+                                        r[b3] = r[b3]/coeff
+                                        r[a3] = r[a3]/coeff
+                                rn = r[c3].match(a4*a4*x**2-b4*b4)
+                                if check==0: # if r[c3] becomes zero at x0
+                                    rn = r[c3].match(a4*a4*x**2)
+                                    if rn:
+                                        rn[b4] = 0
+                                if rn and r[b3] != 0:
+                                    rn = {'n':rn[b4], 'a4':rn[a4]}
+                                    rn['c4'] = r[b3].match(c4*x)[b4]
+                                    matching_hints["2nd_linear_bessel"] = rn
+
+                #If the ODE is ordinary and is of the form of Airy's Equation
+                #Eq(x**2*Derivative(y(x),x,x)-(ax+b)*y(x))
+
+                if p.is_zero:
+                    a4 = Wild('a4', exclude=[x,f(x),df])
+                    b4 = Wild('b4', exclude=[x,f(x),df])
+                    rn = q.match(a4+b4*x)
+                    if rn and rn[b4] != 0:
+                        rn = {'b':rn[a4],'m':rn[b4]}
+                        matching_hints["2nd_linear_airy"] = rn
     if order > 0:
         # Any ODE that can be solved with a substitution and
         # repeated integration e.g.:
@@ -1445,7 +1478,6 @@ def _test_term(coeff, order):
 
     # Order keys based on allhints.
     retlist = [i for i in allhints if i in matching_hints]
-
     if dict:
         # Dictionaries are ordered arbitrarily, so make note of which
         # hint would come first for dsolve().  Use an ordered dict in Py 3.
@@ -3841,6 +3873,37 @@ def ode_2nd_power_series_ordinary(eq, func, order, match):
     return Eq(f(x), series)
 
 
+def ode_2nd_linear_airy(eq, func, order, match):
+    r"""
+    Gives solution of the Airy differential equation
+
+    .. math :: \frac{d^2y}{dx^2} + (a + b x) y(x) = 0
+
+    in terms of Airy special functions airyai and airybi.
+
+    Examples
+    ========
+
+    >>> from sympy import dsolve, Function, pprint
+    >>> from sympy.abc import x
+    >>> f = Function("f")
+    >>> eq = f(x).diff(x, 2) - x*f(x)
+    >>> dsolve(eq)
+    Eq(f(x), C1*airyai(x) + C2*airybi(x))
+    """
+    x = func.args[0]
+    f = func.func
+    C0, C1 = get_numbered_constants(eq, num=2)
+    b = match['b']
+    m = match['m']
+    if m.is_positive:
+        arg = - b/cbrt(m)**2 - cbrt(m)*x
+    elif m.is_negative:
+        arg = - b/cbrt(-m)**2 + cbrt(-m)*x
+    else:
+        arg = - b/cbrt(-m)**2 + cbrt(-m)*x
+    return Eq(f(x), C0*airyai(arg) + C1*airybi(arg))
+
 
 def ode_2nd_power_series_regular(eq, func, order, match):
     r"""
@@ -3964,6 +4027,40 @@ def ode_2nd_power_series_regular(eq, func, order, match):
             return Eq(f(x), collect(finalseries1 + finalseries2,
                 [C0, C1]) + Order(x**terms))
 
+def ode_2nd_linear_bessel(eq, func, order, match):
+    r"""
+    Gives solution of the Bessel differential equation
+
+    .. math :: x**2*\frac{d^2y}{dx^2} + x*\frac{dy}{dx}*y(x) + (x**2-n**2)*y(x)
+
+    if n is integer then the solution is of the form Eq(f(x), C0*besselj(n,x)
+    + C1*bessely(n,x)) as both the solutions are linearly independant else if
+    n is a fraction then the solution is of the form Eq(f(x), C0*besselj(n,x)
+    + C1*besselj(-n,x)) which can also transform into Eq(f(x), C0*besselj(n,x)
+    + C1*bessely(n,x)).
+
+    >>> from sympy.abc import x, y, a, v
+    >>> from sympy.solvers.ode import dsolve, checkodesol
+    >>> from sympy import pprint, Function
+    >>> f = Function('f')
+    >>> y = f(x)
+    >>> genform = x**2*y.diff(x, 2) + x*y.diff(x) + (x**2 - v**2)*y
+    >>> dsolve(genform)
+    Eq(f(x), C1*besselj(sqrt(v**2), x) + C2*bessely(sqrt(v**2), x))
+
+
+    https://www.math24.net/bessel-differential-equation/
+
+    """
+    x = func.args[0]
+    f = func.func
+    C0, C1 = get_numbered_constants(eq, num=2)
+    n = match['n']
+    a4 = match['a4']
+    c4 = match['c4']
+    n = sqrt(n**2 + Rational(1, 4)*(c4 - 1)**2)
+    return (Eq(f(x), (x**(Rational(1-c4,2)))*(C0*besselj(n,x) + C1*bessely(n,x)).subs(x,a4*x)))
+
 def _frobenius(n, m, p0, q0, p, q, x0, x, c, check=None):
     r"""
     Returns a dict with keys as coefficients and values as their values in terms of C0
diff --git a/sympy/solvers/tests/test_ode.py b/sympy/solvers/tests/test_ode.py
index b5fb00882f08..61f49a2fdba9 100644
--- a/sympy/solvers/tests/test_ode.py
+++ b/sympy/solvers/tests/test_ode.py
@@ -6,6 +6,9 @@
     checkodesol, classify_ode, classify_sysode, constant_renumber,
     constantsimp, homogeneous_order, infinitesimals, checkinfsol,
     checksysodesol, solve_ics, dsolve, get_numbered_constants)
+
+from sympy.functions import airyai, airybi, besselj, bessely
+
 from sympy.solvers.deutils import ode_order
 from sympy.utilities.pytest import XFAIL, skip, raises, slow, ON_TRAVIS
 
@@ -256,14 +259,13 @@ def test_linear_2eq_order2():
     # FIXME: assert checksysodesol(eq7, sol7) == (True, [0, 0])
 
     eq8 = (Eq(diff(x(t),t,t), t*(4*x(t) + 9*y(t))), Eq(diff(y(t),t,t), t*(12*x(t) - 6*y(t))))
-    sol8 = ("[Eq(x(t), -sqrt(133)*((-sqrt(133) - 1)*(C2*(133*t**8/24 - t**3/6 + sqrt(133)*t**3/2 + 1) + "
-    "C1*t*(sqrt(133)*t**4/6 - t**3/12 + 1) + O(t**6)) - (-1 + sqrt(133))*(C2*(-sqrt(133)*t**3/6 - t**3/6 + 1) + "
-    "C1*t*(-sqrt(133)*t**3/12 - t**3/12 + 1) + O(t**6)) - 4*C2*(133*t**8/24 - t**3/6 + sqrt(133)*t**3/2 + 1) + "
-    "4*C2*(-sqrt(133)*t**3/6 - t**3/6 + 1) - 4*C1*t*(sqrt(133)*t**4/6 - t**3/12 + 1) + "
-    "4*C1*t*(-sqrt(133)*t**3/12 - t**3/12 + 1) + O(t**6))/3192), Eq(y(t), -sqrt(133)*(-C2*(133*t**8/24 - t**3/6 + "
-    "sqrt(133)*t**3/2 + 1) + C2*(-sqrt(133)*t**3/6 - t**3/6 + 1) - C1*t*(sqrt(133)*t**4/6 - t**3/12 + 1) + "
-    "C1*t*(-sqrt(133)*t**3/12 - t**3/12 + 1) + O(t**6))/266)]")
-    assert str(dsolve(eq8)) == sol8
+    sol8 = [Eq(x(t), -sqrt(133)*(-4*C1*airyai(t*(-1 + sqrt(133))**(S(1)/3)) + 4*C1*airyai(-t*(1 + \
+    sqrt(133))**(S(1)/3)) - 4*C2*airybi(t*(-1 + sqrt(133))**(S(1)/3)) + 4*C2*airybi(-t*(1 + sqrt(133))**(S(1)/3)) +\
+    (-sqrt(133) - 1)*(C1*airyai(t*(-1 + sqrt(133))**(S(1)/3)) + C2*airybi(t*(-1 + sqrt(133))**(S(1)/3))) - (-1 +\
+    sqrt(133))*(C1*airyai(-t*(1 + sqrt(133))**(S(1)/3)) + C2*airybi(-t*(1 + sqrt(133))**(S(1)/3))))/3192), \
+    Eq(y(t), -sqrt(133)*(-C1*airyai(t*(-1 + sqrt(133))**(S(1)/3)) + C1*airyai(-t*(1 + sqrt(133))**(S(1)/3)) -\
+    C2*airybi(t*(-1 + sqrt(133))**(S(1)/3)) + C2*airybi(-t*(1 + sqrt(133))**(S(1)/3)))/266)]
+    assert dsolve(eq8) == sol8
     # FIXME: assert checksysodesol(eq8, sol8) == (True, [0, 0])
 
     eq9 = (Eq(diff(x(t),t,t), t*(4*diff(x(t),t) + 9*diff(y(t),t))), Eq(diff(y(t),t,t), t*(12*diff(x(t),t) - 6*diff(y(t),t))))
@@ -2496,7 +2498,6 @@ def test_nth_order_linear_euler_eq_nonhomogeneous_undetermined_coefficients():
     assert our_hint in classify_ode(eq, f(x))
     assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
     assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
-
     eq = Eq(x**2*diff(f(x), x, x) + 3*x*diff(f(x), x) - 8*f(x), log(x)**3 - log(x))
     sol = C1/x**4 + C2*x**2 - Rational(1,8)*log(x)**3 - Rational(3,32)*log(x)**2 - Rational(1,64)*log(x) - Rational(7, 256)
     sols = constant_renumber(sol)
@@ -2978,15 +2979,14 @@ def test_2nd_power_series_ordinary():
     # checkodesol doesn't work with them.
     C1, C2 = symbols("C1 C2")
     eq = f(x).diff(x, 2) - x*f(x)
-    assert classify_ode(eq) == ('2nd_power_series_ordinary',)
-    assert dsolve(eq) == Eq(f(x),
+    assert classify_ode(eq) == ('2nd_linear_airy', '2nd_power_series_ordinary')
+    assert dsolve(eq, hint='2nd_power_series_ordinary') == Eq(f(x),
         C2*(x**3/6 + 1) + C1*x*(x**3/12 + 1) + O(x**6))
-    assert dsolve(eq, x0=-2) == Eq(f(x),
+    assert dsolve(eq, x0=-2, hint='2nd_power_series_ordinary') == Eq(f(x),
         C2*((x + 2)**4/6 + (x + 2)**3/6 - (x + 2)**2 + 1)
         + C1*(x + (x + 2)**4/12 - (x + 2)**3/3 + S(2))
         + O(x**6))
-    assert dsolve(eq, n=2) == Eq(f(x), C2*x + C1 + O(x**2))
-
+    assert dsolve(eq, n=2, hint='2nd_power_series_ordinary') == Eq(f(x), C2*x + C1 + O(x**2))
     eq = (1 + x**2)*(f(x).diff(x, 2)) + 2*x*(f(x).diff(x)) -2*f(x)
     assert classify_ode(eq) == ('2nd_power_series_ordinary',)
     assert dsolve(eq) == Eq(f(x), C2*(-x**4/3 + x**2 + 1) + C1*x
@@ -3004,10 +3004,24 @@ def test_2nd_power_series_ordinary():
         + 1) + O(x**6))
 
     eq = f(x).diff(x, 2) + x*f(x)
-    assert classify_ode(eq) == ('2nd_power_series_ordinary',)
-    assert dsolve(eq, n=7) == Eq(f(x), C2*(
+    assert classify_ode(eq) == ('2nd_linear_airy', '2nd_power_series_ordinary')
+    assert dsolve(eq, n=7, hint='2nd_power_series_ordinary') == Eq(f(x), C2*(
         x**6/180 - x**3/6 + 1) + C1*x*(-x**3/12 + 1) + O(x**7))
 
+def test_Airy_equation():
+    C1, C2 = symbols("C1 C2")
+    eq = f(x).diff(x, 2) - x*f(x)
+    sol = dsolve(eq)
+    assert classify_ode(eq) == ("2nd_linear_airy",'2nd_power_series_ordinary')
+    assert sol == Eq(f(x), C1*airyai(x) + C2*airybi(x))
+    assert checkodesol(eq, sol) == (True, 0)
+    eq = f(x).diff(x, 2) + 2*x*f(x)
+    assert classify_ode(eq) == ("2nd_linear_airy",'2nd_power_series_ordinary')
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*airyai(-2**(S(1)/3)*x) + C2*airybi(-2**(S(1)/3)*x))
+    assert checkodesol(eq, sol) == (True, 0)
+
+
 
 def test_2nd_power_series_regular():
     # FIXME: Maybe there should be a way to check series solutions
@@ -3029,12 +3043,35 @@ def test_2nd_power_series_regular():
         + O(x**6))
 
     eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 - S(1)/4)*f(x)
-    assert dsolve(eq) == Eq(f(x), C1*(x**4/24 - x**2/2 + 1)/sqrt(x) +
+    assert dsolve(eq, hint='2nd_power_series_regular') == Eq(f(x), C1*(x**4/24 - x**2/2 + 1)/sqrt(x) +
         C2*sqrt(x)*(x**4/120 - x**2/6 + 1) + O(x**6))
 
     eq = x*(f(x).diff(x, 2)) - f(x).diff(x) + 4*x**3*f(x)
     assert dsolve(eq) == Eq(f(x), C2*(-x**4/2 + 1) + C1*x**2 + O(x**6))
 
+def test_2nd_linear_bessel_equation():
+    C1, C2 = symbols("C1 C2")
+    eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 - 4)*f(x)
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*besselj(2, x) + C2*bessely(2, x))
+    assert checkodesol(eq, sol) == (True, 0)
+
+    eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 +25)*f(x)
+    assert classify_ode(eq) == ('2nd_linear_bessel', '2nd_power_series_regular')
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*besselj(5*I, x) + C2*bessely(5*I, x))
+    #assert checkodesol(eq, sol) == (True, 0)
+
+    eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2)*f(x)
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*besselj(0, x) + C2*bessely(0, x))
+    assert checkodesol(eq, sol) == (True, 0)
+
+    eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (81*x**2 -S(1)/9)*f(x)
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*besselj(S(1)/3, 9*x) + C2*bessely(S(1)/3, 9*x))
+    #assert checkodesol(eq, sol) == (True, 0)
+
 def test_issue_7093():
     x = Symbol("x") # assuming x is real leads to an error
     sol = [Eq(f(x), C1 - 2*x*sqrt(x**3)/5),


Code review comments:
diff_hunk: @@ -307,6 +307,9 @@
     "nth_linear_constant_coeff_variation_of_parameters",
     "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters",
     "Liouville",
+    "order_reducible",
path: sympy/solvers/ode.py
text: @user1:
We have nth order reducible below. This shouldn't be added.

diff_hunk: @@ -3029,12 +3043,35 @@ def test_2nd_power_series_regular():
         + O(x**6))
 
     eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 - S(1)/4)*f(x)
-    assert dsolve(eq) == Eq(f(x), C1*(x**4/24 - x**2/2 + 1)/sqrt(x) +
+    assert dsolve(eq, hint='2nd_power_series_regular') == Eq(f(x), C1*(x**4/24 - x**2/2 + 1)/sqrt(x) +
         C2*sqrt(x)*(x**4/120 - x**2/6 + 1) + O(x**6))
 
     eq = x*(f(x).diff(x, 2)) - f(x).diff(x) + 4*x**3*f(x)
     assert dsolve(eq) == Eq(f(x), C2*(-x**4/2 + 1) + C1*x**2 + O(x**6))
 
+def test_2nd_linear_bessel_equation():
+    C1, C2 = symbols("C1 C2")
+    eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 - 4)*f(x)
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*besselj(2, x) + C2*bessely(2, x))
+    assert checkodesol(eq, sol) == (True, 0)
+
+    eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 +25)*f(x)
+    assert classify_ode(eq) == ('2nd_linear_bessel', '2nd_power_series_regular')
+    sol = dsolve(eq)
+    assert sol == Eq(f(x), C1*besselj(5*I, x) + C2*bessely(5*I, x))
+    #assert checkodesol(eq, sol) == (True, 0)
path: sympy/solvers/tests/test_ode.py
text: @user1:
Does checkodesol not work here?

@author:
No

@user1:
#16688

diff_hunk: @@ -256,14 +259,13 @@ def test_linear_2eq_order2():
     # FIXME: assert checksysodesol(eq7, sol7) == (True, [0, 0])
 
     eq8 = (Eq(diff(x(t),t,t), t*(4*x(t) + 9*y(t))), Eq(diff(y(t),t,t), t*(12*x(t) - 6*y(t))))
-    sol8 = ("[Eq(x(t), -sqrt(133)*((-sqrt(133) - 1)*(C2*(133*t**8/24 - t**3/6 + sqrt(133)*t**3/2 + 1) + "
-    "C1*t*(sqrt(133)*t**4/6 - t**3/12 + 1) + O(t**6)) - (-1 + sqrt(133))*(C2*(-sqrt(133)*t**3/6 - t**3/6 + 1) + "
-    "C1*t*(-sqrt(133)*t**3/12 - t**3/12 + 1) + O(t**6)) - 4*C2*(133*t**8/24 - t**3/6 + sqrt(133)*t**3/2 + 1) + "
-    "4*C2*(-sqrt(133)*t**3/6 - t**3/6 + 1) - 4*C1*t*(sqrt(133)*t**4/6 - t**3/12 + 1) + "
-    "4*C1*t*(-sqrt(133)*t**3/12 - t**3/12 + 1) + O(t**6))/3192), Eq(y(t), -sqrt(133)*(-C2*(133*t**8/24 - t**3/6 + "
-    "sqrt(133)*t**3/2 + 1) + C2*(-sqrt(133)*t**3/6 - t**3/6 + 1) - C1*t*(sqrt(133)*t**4/6 - t**3/12 + 1) + "
-    "C1*t*(-sqrt(133)*t**3/12 - t**3/12 + 1) + O(t**6))/266)]")
-    assert str(dsolve(eq8)) == sol8
+    sol8 = [Eq(x(t), -sqrt(133)*(-4*C1*airyai(t*(-1 + sqrt(133))**(S(1)/3)) + 4*C1*airyai(-t*(1 + \
+    sqrt(133))**(S(1)/3)) - 4*C2*airybi(t*(-1 + sqrt(133))**(S(1)/3)) + 4*C2*airybi(-t*(1 + sqrt(133))**(S(1)/3)) +\
+    (-sqrt(133) - 1)*(C1*airyai(t*(-1 + sqrt(133))**(S(1)/3)) + C2*airybi(t*(-1 + sqrt(133))**(S(1)/3))) - (-1 +\
+    sqrt(133))*(C1*airyai(-t*(1 + sqrt(133))**(S(1)/3)) + C2*airybi(-t*(1 + sqrt(133))**(S(1)/3))))/3192), \
+    Eq(y(t), -sqrt(133)*(-C1*airyai(t*(-1 + sqrt(133))**(S(1)/3)) + C1*airyai(-t*(1 + sqrt(133))**(S(1)/3)) -\
+    C2*airybi(t*(-1 + sqrt(133))**(S(1)/3)) + C2*airybi(-t*(1 + sqrt(133))**(S(1)/3)))/266)]
+    assert dsolve(eq8) == sol8
     # FIXME: assert checksysodesol(eq8, sol8) == (True, [0, 0])
path: sympy/solvers/tests/test_ode.py
text: @user1:
Can you use checksysodesol here?

@author:

> Can you add something like the code from #16688 to checkodesol So that we can use it with all the Bessel function solutions here?



  okk I will try in new pr .


@user1:
Or otherwise I can do it if you prefer

@author:
I have made a pr for this you can check #16707.