DECLARE SUB RK4 (l%, x#, y#, z#, h#, x1#, y1#, z1#)
DECLARE FUNCTION F# (x#, y#, z#)
DECLARE FUNCTION G# (l%, x#, y#, z#)
DECLARE FUNCTION FX# (x#)
'**************************************************************************
'*  Solve a boundary value problem for a second order DE of the form:     *
'*  y" - p(x)*y = f(x) with conditions: y(x0)=a and y(x1)=b (x0<=x<=x1)   *
'*  using a Runge-Kutta integration method.                               *
'* ---------------------------------------------------------------------- *
'* Explanations: if y0(x) is a particular solution of y"-p(x)y = f(x) (1) *
'* and y1(x), y2(x) two lineary independant solutions of y"-p(x)y = 0 (2) *
'* (right hand side removed), then the general solution of (1) has the    *
'* form:  y(x) = y0(x) + c1*y1(x) + c2*y2(x), the constants c1, c2 are    *
'* obtained from boundary conditions.                                     *
'* A more efficient way to proceed is: 1) calculate a particular solution *
'* y1(x) of (2) such as y1(x0)=0 and 2) calculate a particular solution   *
'* y2(x) of (1) such as y2(x0)=a, then the general solution of (1) has the*
'* form:    y(x) = y1(x) + C * y2(x), constant C is determined by:        *
'*          C * y2(x1) + y1(x1) = b                                       *
'* ---------------------------------------------------------------------- *
'* SAMPLE RUN:                                                            *
'* Integrate y" - (1/x)y' - (3/xý)y = 5xý from x0=1 to x1=2 with boundary *
'* conditions: y(x0)=11 and y(x1)=36. The exact solution is:              *
'*           y(x) = 8/x + 2xý + x4                                        *
'* The initial 2nd order problem is replaced by the following first order *
'* DE systems:                                                            *
'*                 1st system                2nd system                   *
'*                 y' = z                    y' = z                       *
'*                 z' = z/xý + 3y/xý         z' = z/xý + 3y/xý + 5xý      *
'*                 with z'(x0)=0             with z'(x0)=11               *
'*                                                                        *
'*     X     Y estimated    Y true      Absolute Error                    *
'* ----------------------------------------------------                   *
'*   1.0000  11.00000000   11.00000000    0.00000000                      *
'*   1.0500  11.14980444   11.14980387    0.00000057                      *
'*   1.1000  11.39882827   11.39882727    0.00000100                      *
'*   1.1500  11.74727930   11.74727799    0.00000131                      *
'*   1.2000  12.19626820   12.19626667    0.00000153                      *
'*   1.2500  12.74765794   12.74765625    0.00000169                      *
'*   1.3000  13.40394794   13.40394615    0.00000178                      *
'*   1.3500  14.16818401   14.16818218    0.00000184                      *
'*   1.4000  15.04388756   15.04388571    0.00000185                      *
'*   1.4500  16.03499946   16.03499763    0.00000183                      *
'*   1.5000  17.14583511   17.14583333    0.00000178                      *
'*   1.5500  18.38104827   18.38104657    0.00000170                      *
'*   1.6000  19.74560159   19.74560000    0.00000159                      *
'*   1.6500  21.24474257   21.24474110    0.00000147                      *
'*   1.7000  22.88398367   22.88398235    0.00000132                      *
'*   1.7500  24.66908597   24.66908482    0.00000115                      *
'*   1.8000  26.60604540   26.60604444    0.00000096                      *
'*   1.8500  28.70108132   28.70108057    0.00000075                      *
'*   1.9000  30.96062683   30.96062632    0.00000052                      *
'*   1.9500  33.39132062   33.39132035    0.00000027                      *
'*   2.0000  36.00000000   36.00000000    0.00000000                      *
'*                                                                        *
'* ---------------------------------------------------------------------- *
'* REFERENCE: "Methode de calcul numerique- Tome 2 - Programmes en Basic  *
'*             et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"*
'*             [BIBLI 04].                                                *
'*                                                                        *
'*                             Quick Basic Release By J-P Moreau, Paris.  *
'*                                         (www.jpmoreau.fr)              *
'**************************************************************************
'Program Limits
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0


  ISIZE = 100

DIM x(ISIZE), y(3, ISIZE), z(3, ISIZE)

  x0 = 1#         'starting x
  xl = 2#         'ending x
  kl = 20         'number of steps in x
  h = (xl - x0) / kl 'integration step
  ya = 11#        'starting value for y
  yb = 36#        'final value for y''
  x(0) = x0
  'integration loop
  FOR l = 1 TO 2
    y(l, 0) = (l - 1) * ya: z(l, 0) = 1#
    FOR k = 0 TO kl - 1
      RK4 l, x(k), y(l, k), z(l, k), h, x(k + 1), y(l, k + 1), z(l, k + 1)
    NEXT k
  NEXT l
  'solution is y0=y2+cc*y1, cc is such as y(x1)=yb
  cc = (yb - y(2, kl)) / y(1, kl)
  'write header
  CLS
  PRINT
  PRINT "     X     Y estimated    Y true      Absolute Error "
  PRINT " ----------------------------------------------------"
  'write kl+1 result lines
  FF$ = "  ###.########"
  FOR k = 0 TO kl
    y(0, k) = y(2, k) + cc * y(1, k)
    yex = FX(x(k))        'exact value
    er = ABS(yex - y(0, k))'current error
    PRINT USING "    #.##"; x(k);
    PRINT USING FF$; y(0, k);
    PRINT USING FF$; yex;
    PRINT USING FF$; er
  NEXT k
  PRINT " ----------------------------------------------------"

END 'of main program


  'y'=z
  FUNCTION F (x, y, z)
    F = z
  END FUNCTION

  'exact solution: y=8/x + 2x3 + x4
  FUNCTION FX (x)
    FX = 8# / x + x * x * x * (x + 2#)
  END FUNCTION

  'l=1: z'=z/x²+3y/x² or l=2: z'=z/x²+3y/x²+5x²
  FUNCTION G (l, x, y, z)
    gg = z / x + 3! * y / (x * x)
    IF l = 2 THEN gg = gg + 5# * x * x
    G = gg
  END FUNCTION

  'Integrate sytem from x to x+h using Runge-Kutta
  SUB RK4 (l, x, y, z, h, x1, y1, z1)
    c1 = F(x, y, z)
    d1 = G(l, x, y, z)
    h2 = h / 2#
    c2 = F(x + h2, y + h2 * c1, z + h2 * d1)
    d2 = G(l, x + h2, y + h2 * c1, z + h2 * d1)
    c3 = F(x + h2, y + h2 * c2, z + h2 * d2)
    d3 = G(l, x + h2, y + h2 * c2, z + h2 * d2)
    c4 = F(x + h, y + h * c3, z + h * d3)
    d4 = G(l, x + h, y + h * c3, z + h * d3)
    x1 = x + h
    y1 = y + h * (c1 + 2# * c2 + 2# * c3 + c4) / 6#
    z1 = z + h * (d1 + 2# * d2 + 2# * d3 + d4) / 6#
  END SUB