DECLARE SUB RK4n (n%, x#, y#(), h#, y1#())
DECLARE SUB F (n%, x#, y#(), yp#())
'**************************************************************************
'* Solve an ordinary system of first order differential equations (N<=10) *
'* with initial conditions using a Runge-Kutta integration method.        *
'* ---------------------------------------------------------------------- *
'* SAMPLE RUN:                                                            *
'*  Integrate first order DE system:                                      *
'*  y1' = y2                                                              *
'*  y2' = y3                                                              *
'*  y3' = y4                                                              *
'*  y4' = y5                                                              *
'*  y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)         *
'*  from x=0 to x=2 with initial conditions: y(0)..y(4) = 1               *
'*  (25 integration steps).                                               *
'*                                                                        *
'*     X         Y1 <------------------------> YN estimated               *
'* ---------------------------------------------------------------------- *
'*   0.0000   1.0000000   1.0000000   1.0000000   1.0000000   1.0000000   *
'*   0.0800   1.0832871   1.0832862   1.0832443   1.0816193   1.0382909   *
'*   0.1600   1.1735104   1.1734962   1.1731252   1.1657018   1.0606156   *
'*   0.2400   1.2712454   1.2711671   1.2697821   1.2507235   1.0606138   *
'*   0.3200   1.3771106   1.3768421   1.3732111   1.3346056   1.0307873   *
'*   0.4000   1.4917679   1.4910565   1.4832168   1.4146287   0.9626611   *
'*   0.4800   1.6159210   1.6143206   1.5993581   1.4873681   0.8471623   *
'*   0.5600   1.7503129   1.7470979   1.7208913   1.5486703   0.6752862   *
'*   0.6400   1.8957208   1.8897779   1.8467119   1.5936970   0.4391094   *
'*   0.7200   2.0529492   2.0426461   1.9753055   1.6170637   0.1331626   *
'*   0.8000   2.2228199   2.2058489   2.1047132   1.6131014  -0.2438996   *
'*   0.8800   2.4061600   2.3793580   2.2325243   1.5762584  -0.6875083   *
'*   0.9600   2.6037859   2.5629350   2.3559061   1.5016398  -1.1855881   *
'*   1.0400   2.8164856   2.7561004   2.4716817   1.3856519  -1.7176006   *
'*   1.1200   3.0449981   2.9581109   2.5764618   1.2266816  -2.2547153   *
'*   1.2000   3.2899923   3.1679497   2.6668279   1.0257030  -2.7615081   *
'*   1.2800   3.5520448   3.3843326   2.7395571   0.7866810  -3.1993222   *
'*   1.3600   3.8316195   3.6057337   2.7918674   0.5166480  -3.5309960   *
'*   1.4400   4.1290502   3.8304300   2.8216525   0.2253716  -3.7261918   *
'*   1.5200   4.4445267   4.0565635   2.8276730  -0.0753912  -3.7662105   *
'*   1.6000   4.7780875   4.2822163   2.8096752  -0.3729657  -3.6471428   *
'*   1.6800   5.1296178   4.5051904   2.7684146  -0.6549877  -3.3805446   *
'*   1.7600   5.4988544   4.7245867   2.7055842  -0.9105780  -2.9914571   *
'*   1.8400   5.8853965   4.9378740   2.6236555  -1.1312684  -2.5142924   *
'*   1.9200   6.2887215   5.1439429   2.5256617  -1.3115454  -1.9876112   *
'*   2.0000   6.7082039   5.3416408   2.4149534  -1.4489720  -1.4489720   *
'*                                                                        *
'* ---------------------------------------------------------------------- *
'* REFERENCE: "Methode de calcul numerique- Tome 2 - Programmes en Basic  *
'*             et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"*
'*             [BIBLI 04].                                                *
'*                                                                        *
'*                                  Basic Release By J-P Moreau, Paris.   *
'*                                           (www.jpmoreau.fr)            *
'**************************************************************************
'PROGRAM Test_Rk4n
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

  ISIZE = 9        'maximum number of DEs
 
  DIM y(ISIZE), y1(ISIZE)
  DIM finesse AS INTEGER

  n = 5          ' size of DE system
  x0 = 0#        ' starting x
  xl = 2#        ' ending x
  kl = 25        ' number of tabulated values in x
  finesse = 20
  h = (xl - x0) / kl / finesse' elementary integration step
  x = x0
  FOR i = 0 TO n - 1
    y(i) = 1# ' starting values
  NEXT i
  ' write header
  CLS
  PRINT
  PRINT "     X         Y1 <------------------------> YN estimated"
  PRINT " ----------------------------------------------------------------------"
  ' integration loop
  RK4n n, x, y(), h, y1()
  PRINT USING "####.####"; x;
  FOR i = 0 TO n - 1
    PRINT USING "####.#######"; y(i);
  NEXT i
  PRINT
  FOR k = 1 TO kl * finesse
    x = x + h
    FOR i = 0 TO n - 1
      y(i) = y1(i)
    NEXT i
    RK4n n, x, y(), h, y1()
    IF (k MOD finesse) = 0 THEN
      PRINT USING "####.####"; x;
      FOR i = 0 TO n - 1
        PRINT USING "####.#######"; y(i);
      NEXT i
      PRINT
    END IF
  NEXT k
  PRINT " ----------------------------------------------------------------------"
  PRINT

END

' DE system to integrate (n=5)
SUB F (n, x, y(), yp())
  yp(0) = y(1)
  yp(1) = y(2)
  yp(2) = y(3)
  yp(3) = y(4)
  yp(4) = (45# * y(2) * y(3) * y(4) - 40# * y(3) * y(3) * y(3)) / (9# * y(2) * y(2))
END SUB

'Integrate first order DE sytem from x to x+h using Runge-Kutta
SUB RK4n (n, x, y(), h, y1())
  DIM c1(n - 1), c2(n - 1), c3(n - 1), c4(n - 1), yy(n - 1)
  F n, x, y(), c1()
  h2 = h / 2#
  FOR i = 0 TO n - 1
    yy(i) = y(i) + h2 * c1(i)
  NEXT i
  F n, x + h2, yy(), c2()
  FOR i = 0 TO n - 1
    yy(i) = y(i) + h2 * c2(i)
  NEXT i
  F n, x + h2, yy(), c3()
  FOR i = 0 TO n - 1
    yy(i) = y(i) + h * c3(i)
  NEXT i
  F n, x + h, yy(), c4()
  FOR i = 0 TO n - 1
    y1(i) = y(i) + h * (c1(i) + 2# * c2(i) + 2# * c3(i) + c4(i)) / 6#
  NEXT i

END SUB