Airy equation

We found the Taylor series of the solution of the initial value problem. Let's write out a reasonably large partial sum of the terms and plot it.

import sympy as sp
sp.init_printing()
x,y0,v0 = sp.symbols('x y0 v0')

def approxe(N):
    term1 = y0  ; sum1 = term1;
    term2 = v0*x; sum2 = term2
    for n in range(1,N):
        term1 *= x**3/(3*n  )/(3*n-1); sum1 += term1
        term2 *= x**3/(3*n+1)/(3*n  ); sum2 += term2
    return sum1 + sum2

def approx(N,Y0,V0):
    return sp.lambdify(x,approxe(N).subs({y0:Y0,v0:V0}),'numpy')

approxe(15)

$$\frac{v_{0} x^{43}}{476620138022142319989405253632000000} + \frac{v_{0} x^{40}}{263909268007830741965340672000000} + \frac{v_{0} x^{37}}{169172607697327398695731200000} + \frac{v_{0} x^{34}}{127006462235230779801600000} + \frac{v_{0} x^{31}}{113196490405731532800000} + \frac{v_{0} x^{28}}{121716656350248960000} + \frac{v_{0} x^{25}}{161000868188160000} + \frac{v_{0} x^{22}}{268334780313600} + \frac{v_{0} x^{19}}{580811212800} + \frac{v_{0} x^{16}}{1698278400} + \frac{v_{0} x^{13}}{7076160} + \frac{v_{0} x^{10}}{45360} + \frac{v_{0} x^{7}}{504} + \frac{v_{0} x^{4}}{12} + v_{0} x + \frac{x^{42} y_{0}}{52854284846177190343870827724800000} + \frac{x^{39} y_{0}}{30693545206839251070772838400000} + \frac{x^{36} y_{0}}{20710894201645918401331200000} + \frac{x^{33} y_{0}}{16437217620353903493120000} + \frac{x^{30} y_{0}}{15565546988971499520000} + \frac{x^{27} y_{0}}{17891433320656896000} + \frac{x^{24} y_{0}}{25486372251648000} + \frac{x^{21} y_{0}}{46170964224000} + \frac{x^{18} y_{0}}{109930867200} + \frac{x^{15} y_{0}}{359251200} + \frac{x^{12} y_{0}}{1710720} + \frac{x^{9} y_{0}}{12960} + \frac{x^{6} y_{0}}{180} + \frac{x^{3} y_{0}}{6} + y_{0}$$

%matplotlib inline
from pylab import plot,ylim,legend,axhline,axvline
from numpy import linspace
X = linspace(-8,3,500)
axhline(y=0, color='k', alpha=0.3)
axvline(x=0, color='k', alpha=0.3)
plot(X,approx(15,0,1)(X),label='15 terms')
plot(X,approx(25,0,1)(X),label='25 terms')
ylim(-5,5)
legend(loc='lower right');

axhline(y=0, color='k', alpha=0.3)
axvline(x=0, color='k', alpha=0.3)
plot(X,approx(15,1,0)(X),label='15 terms')
plot(X,approx(25,1,0)(X),label='25 terms')
ylim(-5,5)
legend(loc='lower right');

X = linspace(-15,3,500)
axhline(y=0, color='k', alpha=0.3)
axvline(x=0, color='k', alpha=0.3)
plot(X,approx(25,1,0)(X),label='25 terms')
plot(X,approx(35,1,0)(X),label='35 terms')
plot(X,approx(55,1,0)(X),label='55 terms')
ylim(-1,2)
legend(loc='lower right');