Solution of the differential equation

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Cephes Math Library Release 2.8:  June, 2000
'Copyright by Stephen L. Moshier
'
'Contributors:
'    * Sergey Bochkanov (ALGLIB project). Translation from C to
'      pseudocode.
'
'See subroutines comments for additional copyrights.
'
'>>> SOURCE LICENSE >>>
'This program is free software; you can redistribute it and/or modify
'it under the terms of the GNU General Public License as published by
'the Free Software Foundation (www.fsf.org); either version 2 of the
'License, or (at your option) any later version.
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'http://www.fsf.org/licensing/licenses
'
'>>> END OF LICENSE >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Airy function
'
'Solution of the differential equation
'
'y"(x) = xy.
'
'The function returns the two independent solutions Ai, Bi
'and their first derivatives Ai'(x), Bi'(x).
'
'Evaluation is by power series summation for small x,
'by rational minimax approximations for large x.
'
'
'
'ACCURACY:
'Error criterion is absolute when function <= 1, relative
'when function > 1, except * denotes relative error criterion.
'For large negative x, the absolute error increases as x^1.5.
'For large positive x, the relative error increases as x^1.5.
'
'Arithmetic  domain   function  # trials      peak         rms
'IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
'IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
'IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
'IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
'IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
'IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
'
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Sub Airy(ByVal x As Double, _
         ByRef Ai As Double, _
         ByRef Aip As Double, _
         ByRef Bi As Double, _
         ByRef Bip As Double)
    Dim z As Double
    Dim zz As Double
    Dim t As Double
    Dim f As Double
    Dim g As Double
    Dim uf As Double
    Dim ug As Double
    Dim k As Double
    Dim zeta As Double
    Dim theta As Double
    Dim domflg As Long
    Dim c1 As Double
    Dim c2 As Double
    Dim sqrt3 As Double
    Dim sqpii As Double
    Dim AFN As Double
    Dim AFD As Double
    Dim AGN As Double
    Dim AGD As Double
    Dim APFN As Double
    Dim APFD As Double
    Dim APGN As Double
    Dim APGD As Double
    Dim AN As Double
    Dim AD As Double
    Dim APN As Double
    Dim APD As Double
    Dim BN16 As Double
    Dim BD16 As Double
    Dim BPPN As Double
    Dim BPPD As Double

    sqpii = 0.564189583547756
    c1 = 0.355028053887817
    c2 = 0.258819403792807
    sqrt3 = 1.73205080756888
    domflg = 0#
    If x > 25.77 Then
        Ai = 0#
        Aip = 0#
        Bi = MaxRealNumber
        Bip = MaxRealNumber
        Exit Sub
    End If
    If x < -2.09 Then
        domflg = 15#
        t = Sqr(-x)
        zeta = -(2# * x * t / 3#)
        t = Sqr(t)
        k = sqpii / t
        z = 1# / zeta
        zz = z * z
        AFN = -0.131696323418332
        AFN = AFN * zz - 0.626456544431912
        AFN = AFN * zz - 0.693158036036933
        AFN = AFN * zz - 0.279779981545119
        AFN = AFN * zz - 0.04919001326095
        AFN = AFN * zz - 4.06265923594885E-03
        AFN = AFN * zz - 1.59276496239262E-04
        AFN = AFN * zz - 2.77649108155233E-06
        AFN = AFN * zz - 1.67787698489115E-08
        AFD = 1#
        AFD = AFD * zz + 13.3560420706553
        AFD = AFD * zz + 32.6825032795225
        AFD = AFD * zz + 26.73670409415
        AFD = AFD * zz + 9.1870740290726
        AFD = AFD * zz + 1.47529146771666
        AFD = AFD * zz + 0.115687173795188
        AFD = AFD * zz + 4.40291641615211E-03
        AFD = AFD * zz + 7.54720348287414E-05
        AFD = AFD * zz + 4.5185009297058E-07
        uf = 1# + zz * AFN / AFD
        AGN = 1.97339932091686E-02
        AGN = AGN * zz + 0.391103029615688
        AGN = AGN * zz + 1.06579897599596
        AGN = AGN * zz + 0.93916922981665
        AGN = AGN * zz + 0.351465656105548
        AGN = AGN * zz + 6.33888919628925E-02
        AGN = AGN * zz + 5.85804113048388E-03
        AGN = AGN * zz + 2.82851600836737E-04
        AGN = AGN * zz + 6.98793669997261E-06
        AGN = AGN * zz + 8.11789239554389E-08
        AGN = AGN * zz + 3.41551784765924E-10
        AGD = 1#
        AGD = AGD * zz + 9.30892908077442
        AGD = AGD * zz + 19.8352928718312
        AGD = AGD * zz + 15.5646628932865
        AGD = AGD * zz + 5.47686069422975
        AGD = AGD * zz + 0.954293611618962
        AGD = AGD * zz + 8.64580826352392E-02
        AGD = AGD * zz + 4.12656523824223E-03
        AGD = AGD * zz + 1.01259085116509E-04
        AGD = AGD * zz + 1.17166733214414E-06
        AGD = AGD * zz + 4.9183457006293E-09
        ug = z * AGN / AGD
        theta = zeta + 0.25 * PI()
        f = Sin(theta)
        g = Cos(theta)
        Ai = k * (f * uf - g * ug)
        Bi = k * (g * uf + f * ug)
        APFN = 0.185365624022536
        APFN = APFN * zz + 0.886712188052584
        APFN = APFN * zz + 0.987391981747399
        APFN = APFN * zz + 0.401241082318004
        APFN = APFN * zz + 7.10304926289631E-02
        APFN = APFN * zz + 5.90618657995662E-03
        APFN = APFN * zz + 2.33051409401777E-04
        APFN = APFN * zz + 4.08718778289035E-06
        APFN = APFN * zz + 2.48379932900442E-08
        APFD = 1#
        APFD = APFD * zz + 14.7345854687503
        APFD = APFD * zz + 37.542393343549
        APFD = APFD * zz + 31.4657751203046
        APFD = APFD * zz + 10.9969125207299
        APFD = APFD * zz + 1.78885054766999
        APFD = APFD * zz + 0.141733275753663
        APFD = APFD * zz + 5.44066067017226E-03
        APFD = APFD * zz + 9.39421290654511E-05
        APFD = APFD * zz + 5.65978713036027E-07
        uf = 1# + zz * APFN / APFD
        APGN = -3.55615429033082E-02
        APGN = APGN * zz - 0.637311518129436
        APGN = APGN * zz - 1.70856738884312
        APGN = APGN * zz - 1.50221872117317
        APGN = APGN * zz - 0.563606665822103
        APGN = APGN * zz - 0.102101031120217
        APGN = APGN * zz - 9.48396695961445E-03
        APGN = APGN * zz - 4.60325307486781E-04
        APGN = APGN * zz - 1.14300836484517E-05
        APGN = APGN * zz - 1.33415518685547E-07
        APGN = APGN * zz - 5.63803833958894E-10
        APGD = 1#
        APGD = APGD * zz + 9.8586580169613
        APGD = APGD * zz + 21.6401867356586
        APGD = APGD * zz + 17.3130776389749
        APGD = APGD * zz + 6.17872175280829
        APGD = APGD * zz + 1.08848694396321
        APGD = APGD * zz + 9.95005543440888E-02
        APGD = APGD * zz + 4.78468199683887E-03
        APGD = APGD * zz + 1.18159633322839E-04
        APGD = APGD * zz + 1.37480673554219E-06
        APGD = APGD * zz + 5.79912514929148E-09
        ug = z * APGN / APGD
        k = sqpii * t
        Aip = -(k * (g * uf + f * ug))
        Bip = k * (f * uf - g * ug)
        Exit Sub
    End If
    If x >= 2.09 Then
        domflg = 5#
        t = Sqr(x)
        zeta = 2# * x * t / 3#
        g = Exp(zeta)
        t = Sqr(t)
        k = 2# * t * g
        z = 1# / zeta
        AN = 0.346538101525629
        AN = AN * z + 12.0075952739646
        AN = AN * z + 76.2796053615235
        AN = AN * z + 168.089224934631
        AN = AN * z + 159.756391350164
        AN = AN * z + 70.5360906840444
        AN = AN * z + 14.026469116339
        AN = AN * z + 1#
        AD = 0.56759453263877
        AD = AD * z + 14.7562562584847
        AD = AD * z + 84.5138970141475
        AD = AD * z + 177.3180881454
        AD = AD * z + 164.23469287153
        AD = AD * z + 71.4778400825576
        AD = AD * z + 14.0959135607834
        AD = AD * z + 1#
        f = AN / AD
        Ai = sqpii * f / k
        k = -(0.5 * sqpii * t / g)
        APN = 0.613759184814036
        APN = APN * z + 14.7454670787755
        APN = APN * z + 82.0584123476061
        APN = APN * z + 171.184781360976
        APN = APN * z + 159.317847137142
        APN = APN * z + 69.9778599330103
        APN = APN * z + 13.9470856980482
        APN = APN * z + 1#
        APD = 0.334203677749737
        APD = APD * z + 11.1810297306158
        APD = APD * z + 71.172735214786
        APD = APD * z + 158.778084372838
        APD = APD * z + 153.206427475809
        APD = APD * z + 68.675230459278
        APD = APD * z + 13.8498634758259
        APD = APD * z + 1#
        f = APN / APD
        Aip = f * k
        If x > 8.3203353 Then
            BN16 = -0.253240795869364
            BN16 = BN16 * z + 0.575285167332467
            BN16 = BN16 * z - 0.329907036873225
            BN16 = BN16 * z + 0.06444040689482
            BN16 = BN16 * z - 3.82519546641337E-03
            BD16 = 1#
            BD16 = BD16 * z - 7.15685095054035
            BD16 = BD16 * z + 10.6039580715665
            BD16 = BD16 * z - 5.23246636471251
            BD16 = BD16 * z + 0.957395864378384
            BD16 = BD16 * z - 0.055082814716355
            f = z * BN16 / BD16
            k = sqpii * g
            Bi = k * (1# + f) / t
            BPPN = 0.465461162774652
            BPPN = BPPN * z - 1.08992173800494
            BPPN = BPPN * z + 0.638800117371828
            BPPN = BPPN * z - 0.126844349553103
            BPPN = BPPN * z + 7.6248784434211E-03
            BPPD = 1#
            BPPD = BPPD * z - 8.70622787633159
            BPPD = BPPD * z + 13.8993162704553
            BPPD = BPPD * z - 7.14116144616431
            BPPD = BPPD * z + 1.34008595960681
            BPPD = BPPD * z - 7.84273211323342E-02
            f = z * BPPN / BPPD
            Bip = k * t * (1# + f)
            Exit Sub
        End If
    End If
    f = 1#
    g = x
    t = 1#
    uf = 1#
    ug = x
    k = 1#
    z = x * x * x
    Do While t > MachineEpsilon
        uf = uf * z
        k = k + 1#
        uf = uf / k
        ug = ug * z
        k = k + 1#
        ug = ug / k
        uf = uf / k
        f = f + uf
        k = k + 1#
        ug = ug / k
        g = g + ug
        t = Abs(uf / f)
    Loop
    uf = c1 * f
    ug = c2 * g
    If domflg Mod 2# = 0# Then
        Ai = uf - ug
    End If
    If domflg \ 2# Mod 2# = 0# Then
        Bi = sqrt3 * (uf + ug)
    End If
    k = 4#
    uf = x * x / 2#
    ug = z / 3#
    f = uf
    g = 1# + ug
    uf = uf / 3#
    t = 1#
    Do While t > MachineEpsilon
        uf = uf * z
        ug = ug / k
        k = k + 1#
        ug = ug * z
        uf = uf / k
        f = f + uf
        k = k + 1#
        ug = ug / k
        uf = uf / k
        g = g + ug
        k = k + 1#
        t = Abs(ug / g)
    Loop
    uf = c1 * f
    ug = c2 * g
    If domflg \ 4# Mod 2# = 0# Then
        Aip = uf - ug
    End If
    If domflg \ 8# Mod 2# = 0# Then
        Bip = sqrt3 * (uf + ug)
    End If
End Sub