DECLARE SUB RK4 (X#, Y#, Z#, h#, x1#, y1#, z1#)
DECLARE FUNCTION F# (X#, Y#, Z#)
DECLARE FUNCTION G# (X#, Y#, Z#)
'**************************************************************************
'*  Solve a first order DE system (N=2) of the form:                      *
'*  y' = F(x,y,z), z'=G(x,y,z) with conditions: y(x0)=a and z(x0)=b from  *
'*  x=x0 to x=x1 using a Runge-Kutta integration method.                  *
'* ---------------------------------------------------------------------- *
'* SAMPLE RUN:                                                            *
'*      Integrate first order DE systemy:                                 *
'*          y' = y*z + cos(x) - 0.5*sin(2*x)                              *
'*          z' = y*y + z*z - (1 + sin(x))                                 *
'*      from x=0 to x=2 with initial conditions: y(0)=0 and z(0)=0        *
'*           (25 integration steps).                                      *
'*                                                                        *
'*     X     Y estimated  Z estimated                                     *
'* ------------------------------------                                   *
'*   0.0000   0.0000000    0.0000000                                      *
'*   0.0800   0.0765518   -0.0828580                                      *
'*   0.1600   0.1452904   -0.1700792                                      *
'*   0.2400   0.2050556   -0.2597237                                      *
'*   0.3200   0.2550329   -0.3500533                                      *
'*   0.4000   0.2948004   -0.4396008                                      *
'*   0.4800   0.3243379   -0.5271915                                      *
'*   0.5600   0.3440043   -0.6119278                                      *
'*   0.6400   0.3544874   -0.6931461                                      *
'*   0.7200   0.3567373   -0.7703644                                      *
'*   0.8000   0.3518908   -0.8432315                                      *
'*   0.8800   0.3411943   -0.9114874                                      *
'*   0.9600   0.3259299   -0.9749377                                      *
'*   1.0400   0.3073506   -1.0334409                                      *
'*   1.1200   0.2866242   -1.0869048                                      *
'*   1.2000   0.2647908   -1.1352886                                      *
'*   1.2800   0.2427308   -1.1786041                                      *
'*   1.3600   0.2211465   -1.2169165                                      *
'*   1.4400   0.2005537   -1.2503410                                      *
'*   1.5200   0.1812830   -1.2790362                                      *
'*   1.6000   0.1634888   -1.3031952                                      *
'*   1.6800   0.1471644   -1.3230351                                      *
'*   1.7600   0.1321612   -1.3387869                                      *
'*   1.8400   0.1182104   -1.3506867                                      *
'*   1.9200   0.1049465   -1.3589674                                      *
'*   2.0000   0.0919314   -1.3638530                                      *
'*                                                                        *
'* ---------------------------------------------------------------------- *
'* 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_Rk4
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

  ISIZE = 100

DIM X(ISIZE), Y(ISIZE), Z(ISIZE)

  CLS
  x0 = 0#              'starting x
  xl = 2#              'ending x
  kl = 25              'number of steps in x
  h = (xl - x0) / kl   'integration step
  X(0) = x0
  Y(0) = 0#: Z(0) = 0# 'starting values
  'integration loop
  FOR k = 0 TO kl - 1
    RK4 X(k), Y(k), Z(k), h, X(k + 1), Y(k + 1), Z(k + 1)
  NEXT k
  'write header
  PRINT
  PRINT "     X     Y estimated    Z estimated "
  PRINT " -------------------------------------"
  'write kl+1 result lines
  FF$ = "   ###.#######"
  FOR k = 0 TO kl
    PRINT USING "    #.##"; X(k);
    PRINT USING FF$; Y(k);
    PRINT USING FF$; Z(k)
  NEXT k
  PRINT " -------------------------------------"
  PRINT

END 'of main program

'y'=yz + cos(x) - 0.5sin(2x)
FUNCTION F (X, Y, Z)
  F = Y * Z + COS(X) - .5# * SIN(2# * X)
END FUNCTION

'z'=yy + zz -(1+sin(x))
FUNCTION G (X, Y, Z)
  G = Y * Y + Z * Z - (1# + SIN(X))
END FUNCTION

'Integrate sytem from x to x+h using Runge-Kutta
SUB RK4 (X, Y, Z, h, x1, y1, z1)
  c1 = F(X, Y, Z)
  d1 = G(X, Y, Z)
  h2 = h / 2#
  c2 = F(X + h2, Y + h2 * c1, Z + h2 * d1)
  d2 = G(X + h2, Y + h2 * c1, Z + h2 * d1)
  c3 = F(X + h2, Y + h2 * c2, Z + h2 * d2)
  d3 = G(X + h2, Y + h2 * c2, Z + h2 * d2)
  c4 = F(X + h, Y + h * c3, Z + h * d3)
  d4 = G(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