Differential Equations Manual


 Overview
 
 
 
 
XROM  Function  Desciption
 15,00
 15,01
 15,02
 15,03
 15,04
 15,05
 15,06
 15,07
 15,08
 15,09
 15,10
 15,11
 15,12
 15,13
 15,14
 15,15
 15,16
 15,17
 15,18
 15,19
 15,20
 15,21
 15,22
 15,23
 15,24
 15,25
 15,26
 15,27
 15,28
 15,29
 -DIFFEQTN
  "C41"
  "C51"
  "RK4"
  "RK4B"
  "RK4C"
  "RK4N"
  "RK6"
  "RK6B"
  "RK6N"
  "RK8"
  "RK8B"
  "RK8N"
  "RKF"
  "RKFB"
  "RKFN"
  "RKS2"
  "RKS3"
  "RKSN"
  "RKV"
  "RKVB"
-BULIRSHSTR
  "BST"
  "BST2"
  "BSTN"
  "STRM"
-NTH-ORDER
  "NRK4"
  "SRK4"
  "TRK4"
 Section Header
 Constants for "RK8"
 Constants for "RKV"
 Runge-Kutta 4th order formulae
 RK4 for Systems of 2 Differential Equations
 RK4 for Systems of 3 Differential Equations
 RK4 for Systems of n Differential Equations
 Runge-Kutta 6th order formulae
 RK6 for Systems of 2 Differential Equations
 RK6 for Systems of n Differential Equations
 Runge-Kutta 8th order formulae
 RK8 for Systems of 2 Differential Equations
 RK8 for Systems of n Differential Equations
 Runge-Kutta-Felhberg 4-5 method
 RKF for Systems of 2 Differential Equations
 RKF for Systems of n Differential Equations
 RK4 for Systems of 2 Second order Conservative Eq
 RK4 for Systems of 3 Second order Conservative Eq
 RK4 for Systems of n Second order Conservative Eq
 Runge-Kutta-Verner 5-6 method
 RKV for Systems of 2 Differential Equations
 Section Header
 Bulirsh-Stoer Method
 BST  for Systems of 2 Differential Equations
 BST  for Systems of n Differential Equations
 Stoermer's rule for Conservative 2nd-order Equation
 Section Header
 N-th order differential equations ( 4th order method )
 2-nd order differential equations ( 4th order method )
 3-rd order differential equations ( 4th order method )

 

Constants for "RK8"
 

 "C41" stores these 41 constants into registers R11 thru R51
 They are used by an 11-stage method of orde 8 ( RK8 )
 
 
 
  R11 = 0.8273268354 
  R12 = 0.
  R13 = 0.1726731646
  R14 = 9 
  R15= 14
  R16 = 0.25   
  R17 = 0.8961811934
  R18 = -0.2117115009
  R19 = 7
  R20 = 0.06514850915
  R21 = 0.5766714727
  R22 = 0.1855068535 
  R23 = 0.3865246267
  R24 = -0.4634553896
  R25 = 0.3772937693  
  R26 = 0.1996369936
  R27 = 0.09790045616
  R28 = 0.3285172131
  R29 = -0.3497052863
  R30 = -0.03302551131
  R31 = 0.1289862930 
  R32 = -0.01186868389
  R33 = 0.002002165993
  R34 = 0.9576053440
  R35 = -0.6325461607
  R36 = 0.1527777778
  R37 = -0.009086961101
  R38 = 0.03125  
  R39 = 1.062498447
  R40 = -1.810863083
  R41 = 2.031083139
  R42 = -0.6379313502
  R43 = 0.9401094520
  R44 = -2.229158210
  R45 = 7.553840442
  R46 = -7.164951553
  R47 = 2.451380432
  R48 = -0.5512205631
  R49 = 0.05
  R50 = 0.2722222222
  R51 = 2.8125

 

Constants for "RKV"
 

 "C51" stores 51 constants into R19 thru R69
 These constants are used by "RKV" & "RKVB" which employs an 8-stage method with order 5 with a 6th-order error-estimating formula.
 
 
 
  R19 =  1/18 
  R21 = -1/12
  R24 = -2/81
  R28 =  40/33
  R33 = -369/73
  R39 = -8716/891
  R46 =  3015/256
  
  R22 = 1/4
  R25 = 4/27 
  R29 = -4/11
  R34 = 72/73
  R40 = 656/297
  R47 = -9/4  
  
  
  R26 =  8/81
  R30 = -56/11
  R35 = 5380/219
  R41 = 39520/891
  R48 = -4219/78
  
  
  
  R31 = 54/11
  R36 = -12285/584
  R42 = -416/11 
  R49 = 5985/128
  
  
  
  
  R37 = 2695/1752
  R43 = 52/27 
  R50 = -539/384
  
  
  
  
  
  R44 = 0
  R51 = 0
  
  
  
  
  
  
  R52 = 693/3328 
  R20 = 1/18
  R23 = 1/6
  R27 = 2/9
  R32 = 2/3
  R38 = 1
  R45 = 8/9
  R53 = 1
  R54 = 3/80
  R62 = 33/640 
  R55 = 0  
  R63 = 0
  R56 = 4/25
  R64 = -132/325
  R57 = 243/1120
  R65 = 891/2240
  R58 = 77/160 
  R66 = -33/320
  R59 = 73/700
  R67 = -73/700
  R60 = 0
  R68 = 891/8320
  R61 = 0
  R69 = 2/35

 

Runge-Kutta 4th order formulae
 

-We solve    dy/dx = f(x,y)   with the initial value:   y(x0) = y0

RK4 uses the following formulae:

         k1 = h.f(x;y)
         k2 = h.f(x+h/2;y+k1/2)
         k3 = h.f(x+h/2;y+k2/2)
         k4 = h.f(x+h;y+k3)

  and then,   y(x+h) = y(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6      This formula duplicates the Taylor series up to the term in h4
 

Data Registers:           R00:  f name                   ( Registers R00 to R04 are to be initialized before executing "RK4" )

                                      R01 = x0         R03 = h  = stepsize          R05 & R07: temp
                                      R02 = y0         R04 = N = number of steps
Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

     >>>>>       In other words, this subroutine has to change the stack from

      Y = y
      X = x        into       X' = f(x;y)
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)              ( The solution is y = exp ( -x2 ) )
 
 
 01  LBL "T"
 02  *
 03  ST+ X
 04  CHS
 05  RTN

 
 "T"  ASTO 00

  0     STO 01
  1     STO 02    and if we choose  h = 0.1
0.1    STO 03
10     STO 04    XEQ "RK4"   yields          x = 1                                 ( in 20s )
                                X<>Y                    y(1) = 0.367881                   ( the exact result is 1/e =  0.367879441 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 24 = 16 )
 
 

RK4B:  RK4 for Systems of 2 Differential Equations
 

-"RK4B" solves   dy/dx = f(x,y,z)  ,   dz/dx = g(x,y,z)     with the initial values:    y(x0) = y0  ,   z(x0) = z0

-The preceding formulae are now generalized to a system of 2 differential equations:

                         k1 = h.f(x;y;z)                                             l1 = h.g(x;y;z)
                         k2 = h.f(x+h/2;y+k1/2;z+l1/2)                      l2 = h.g(x+h/2;y+k1/2;z+l1/2)
                         k3 = h.f(x+h/2;y+k2/2;z+l2/2)                      l3 = h.g(x+h/2;y+k2/2;z+l2/2)
                         k4 = h.f(x+h;y+k3;z+l3)                               l4 = h.g(x+h;y+k3;z+l3)

  then,   y(x+h) = y(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6
   and   z(x+h) = z(x) + ( l1 + 2.l2 + 2.l3 + l4 ) / 6           duplicate the Taylor series through the terms in h4
 

Data Registers:                 R00:  f name              ( Registers R00 to R05  are to be initialized before executing "RK4B" )

                               R01 = x0         R04 = h  = stepsize              R06 to R09: temp
                               R02 = y0         R05 = N = number of steps
                               R03 = z0
Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                      assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

  >>>>>    In other words, this subroutine has to change the stack from

  Z = z
  Y = y       to     Y' = f(x;y;z)
  X = x               X' = g(x;y;z)
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h)
           X             /          x0+N.h

 
Example:     y(0) = 1 ; y'(0) = 0 ;  y'' + 2.x.y' + 2.y = 0   Evaluate y(1) and y'(1)            (  The solution is y = exp ( -x2 ) ;  y' = -2.x exp ( -x2 )  )

-This second order equation is equivalent to the system:             y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    where z = y'
 
 
 01  LBL "U"
 02  RCL Z
 03  *
 04  +
 05  ST+ X
 06  CHS
 07  RTN

 
"U"  ASTO 00

 0      STO 01
 1      STO 02
 0      STO 03   and if we choose  h = 0.1
0.1    STO 04
10     STO 05    XEQ "RK4B"   yields       x =  1                                        ( in 32s )
                              RDN                      y(1) =  0.367881                            the exact result is  y(1) = 1/e =  0.367879441
                              RDN                      z(1) = -0.735762                            the exact result is  z(1) = -2/e = -0.735758883
 

-Press R/S to continue with the same h and N.
 

N.B:

-To obtain an error estimate, start over again with a smaller stepsize:
  when h is divided by 2 , errors are approximately divided by 16
 

RK4C:  RK4 for Systems of 3 Differential Equations
 

-"RK4C" solves   dy/dx = f(x,y,z,u)  ,   dz/dx = g(x,y,z,u)  ,  du/dx = h(x,y,z,u)   with the initial values:    y(x0) = y0  ,   z(x0) = z0  ,  u(x0) = u0
 

Data Registers:             R00:  f name       ( Registers R00 to R06  are to be initialized before executing "RK4C" )

                               R01 = x0         R05 = h  =  stepsize              R07 to R11: temp
                               R02 = y0         R06 = N = number of steps
                               R03 = z0
                               R04 = u0
Flags: /
Subroutine:  a program calculating f(x;y;z;u) in Z-register , g(x;y;z;u) in Y-register and h(x;y;z;u) in X-register
                      assuming x is in X-register , y is in Y-register , z is in Z-register and u is in T-register upon entry.

      >>>>>    In other words, this subroutine has to change the stack from

  T = u
  Z = z                Z' = f(x;y;z;u)
  Y = y       to     Y' = g(x;y;z;u)
  X = x               X' = h(x;y;z;u)
 
 
      STACK        INPUTS      OUTPUTS
           T             /       u(x0+N.h)
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h)
           X             /          x0+N.h

 
Example:

     dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1   ( whence N = 10 )
     dz/dx = x.( y + z - u )         z(0) = 1
     du/dx = x.y - z.u                u(0) = 2
 
 
  LBL "T"
  R^
  CHS
  STO L
  R^
  ST* Y
  ST+ L
  CLX
  RCL Z
  R^
  ST+ L
  ST* Y
  X<> Z
  ST+ Y
  ST* Z
  X<> T
  LASTX
  * 
  X<>Y
  RTN

 
 "T"   ASTO 00

  0     STO 01
  1     STO 02  STO 03
  2     STO 04
0.1    STO 05
 10    STO 06    XEQ "RK4C"   >>>>      x   =  1                         ( in 53 seconds )
                                                  RDN     y(1) =  0.258209
                                                  RDN     z(1) =  1.157620
                                                  RDN     u(1) =  0.842179

-Press R/S to continue with the same h and N.
 

RK4N:  RK4 for Systems of n Differential Equations
 

-"RK4N" solves the following system with the "classical" 4th order Runge-Kutta method.

     dy1/dx = f1(x,y1, ... ,yn)
     dy2/dx = f2(x,y1, ... ,yn)
             ..................                         with the initial values:     yi(x0)         i = 1 , ... , n

     dyn/dx = fn(x,y1, ... ,yn)

Data Registers:             R00 = subroutine name        ( Registers R00 thru R03 and R10 thru R10+n are to be initialized before executing "RK4N" )

                                        R01 = n  = number of equations = number of functions            R04 thru R09: temp
                                        R02 = h  = stepsize
                                        R03 = N = number of steps

                                        R10 = x0
                                        R11 = y1(x0)                                                                           R11+n thru R10+4n: temp
                                        R12 =  y2(x0)
                                          ...............         ( initial values )

                                        R10+n = yn(x0)
Flags:  F06-F07
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R11+n, ... , R10+2n  respectively
                        with  x, y1, ... , yn in regiters R10, R11, ... , R10+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:     Let's solve the system:

      dy1/dx = - y1y2y3
      dy2/dx =  x ( y1+y2-y3 )
      dy3/dx =  xy1- y2y3

   with the initial values:   y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2

 x ; y1 ; y2 ; y3  will be in R10 ; R11 ; R12 ; R13 respectively, and we have to write a program to compute the 3 functions
 and store the results in R14 ; R15 ; R16   ( just after the "initial-value registers" ).

-For instance, the following subroutine will do the job - let's call it "T"
 
 
 01  LBL "T"
 02  RCL 11
 03  RCL 12
 04  RCL 13
 05  *
 06  *
 07  CHS
 08  STO 14
 09  RCL 11
 10  RCL 12
 11  +
 12  RCL 13
 13  -
 14  RCL 10
 15  *
 16  STO 15
 17  RCL 10
 18  RCL 11
 19  *
 20  RCL 12
 21  RCL 13
 22  *
 23  -
 24  STO 16
 25  RTN

 
We store the name of the subroutine in R00    "T"  ASTO 00

the numbers of equations in R01                      3    STO 01

then the initial values:

     x     in R10                      0 STO 10
     y1   in R11                      1 STO 11
     y2   in R12                      1 STO 12
     y3   in R13                      2 STO 13

-Suppose we want to know   y1(1) , y2(1) ,  y3(1)  with a step size h = 0.1 and therefore N = 10

    h is stored in R02                 0.1   STO 02
   N is stored in R03                 10    STO 03   and  XEQ "RK4N"

about 2mn24s later, the HP-41 returns:

             x  =  1                        in X and R10
        y1(1) =  0.2582093855  in Y and R11
        y2(1) =  1.157619553    in Z and R12
        y3(1) =  0.8421786516  in T and R13

-To continue with the same h and N, simply press R/S and you'll have y1(2) y2(2) y3(2)

-An estimation of the accuracy may be obtained by doing the calculation again with a smaller step size like h=0.05

 It yields:       y1(1) = 0.2582079989    y2(1) = 1.157623732    y3(1) = 0.8421783424

-Theoretically, the results tend to the exact solution as h tends to zero.
 

Runge-Kutta 6th order formulae
 

   dy/dx = f(x,y)   with the initial value:   y(x0) = y0

-The formula duplicates the Taylor series up to the term in h6
-Methods of this order require 7 stages ( 7 evaluations of the function f within each step)
 

Data Registers:                      R00:  f name                 ( Registers R00 to R04 are to be initialized before executing "RK6" )

                               R01 = x0         R03 =  h = the stepsize                   R05: counter
                               R02 = y0         R04 = N = the number of steps      R06 , ...... , R11 = k1 , ..... , k6

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)
 
 
 01  LBL "T"
 02  *
 03  ST+ X
 04  CHS
 05  RTN

 
 "T"  ASTO 00

  0     STO 01
  1     STO 02    and if we choose  h = 0.1
0.1    STO 03
10     STO 04    XEQ "RK6"   yields          x = 1                       ( in 84s )
                                X<>Y                    y(1) = 0.367879436    ( error =  -5 10-9 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 26 = 64 )
 

RK6B:  RK6 for Systems of 2 Differential Equations
 

   dy/dx = f(x,y,z)  ,   dz/dx = g(x,y,z)     with the initial values:    y(x0) = y0  ,   z(x0) = z0
 

Data Registers:                    R00:  f name              ( Registers R00 to R05  are to be initialized before executing "RK6B" )

                               R01 = x0         R04 = h = stepsize                         R06 to R16: temp
                               R02 = y0         R05 = N = number of steps
                               R03 = z0
Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                      assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:          y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    ;    calculate  y(1) and z(1)
 
 
 01  LBL "U"
 02  RCL Z
 03  *
 04  +
 05  ST+ X
 06  CHS
 07  RTN

 
"U"  ASTO 00

 0      STO 01
 1      STO 02
 0      STO 03    and if we choose  h = 0.1
0.1    STO 04
10     STO 05    XEQ "RK6B"   yields      x =  1                                ( in 2mn30s )
                              RDN                      y(1) =  0.367879433             ( y-error = -9 10-9 )
                              RDN                      z(1) = -0.735758865             ( z-error = 17 10-9 )
 

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 26 = 64 )
 

RK6N:  RK6 for Systems of n Differential Equations
 

     dy1/dx = f1(x,y1, ... ,yn)
     dy2/dx = f2(x,y1, ... ,yn)
             ..................                         with the initial values:     yi(x0)         i = 1 , ... , n

     dyn/dx = fn(x,y1, ... ,yn)
 

Data Registers:             R00 = subroutine name       ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RK6N" )

                                        R01 = n  = number of equations = number of functions        R04 thru R13: temp        R14 thru R19: unused
                                        R02 = h  = stepsize
                                        R03 = N = number of steps

                                        R20 = x0
                                        R21 = y1(x0)                                                                           R21+n thru R20+8n: temp
                                        R22 =  y2(x0)
                                          ...............         ( initial values )

                                        R20+n = yn(x0)
Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+n, ... , R20+2n  respectively
                        with  x, y1, ... , yn in regiters R20, R21, ... , R20+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:

      dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1   ( whence N = 10 )
      dz/dx = x.( y + z - u )         z(0) = 1
       du/dx = x.y - z.u                u(0) = 2
 
 
 01  LBL "T"
 02  RCL 21
 03  RCL 22
 04  RCL 23
 05  *
 06  *
 07  CHS
 08  STO 24
 09  RCL 21
 10  RCL 22
 11  +
 12  RCL 23
 13  -
 14  RCL 20
 15  *
 16  STO 25
 17  LASTX
 18  RCL 21
 19  *
 20  RCL 22
 21  RCL 23
 22  *
 23  -
 24  STO 26
 25  RTN

 
 "T"  ASTO 00

 Initial values:

           0   STO 20
           1   STO 21    STO 22
           2   STO 23

 n = 3      STO 01    ( 3 functions )
 h = 0.1   STO 02
 N = 10   STO 03    XEQ "RK6N"  >>>>     x   = 1                     = R20                      ( in 6mn28s )
                                                       RDN    y(1) = 0.258207889  = R21                                                                         y(1) = 0.258207906
                                                       RDN    z(1) = 1.157623948  = R22         the exact results, rounded to 9D are          z(1) = 1.157623981
                                                       RDN    u(1) = 0.842178329  = R23                                                                         u(1) = 0.842178312

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if needed )
 

Runge-Kutta 8th order formulae
 

-Methods of this order require 11 stages. The following one was discovered in 1972 by Cooper and Verner.
-The formula duplicates the Taylor series up to the term in h8
-"RK8" uses 41 coefficients which are of the form ( a + b 211/2 )/c  where a , b , c are integers.

-They are stored in registers R11 thru R51 by "C41"
 

Data Registers:             R00:  f name            ( Registers R00 to R04 are to be initialized before executing "RK8" )

                               R01 = x0         R03 =  h = the stepsize                      R05  to R10 and R52: temp
                               R02 = y0         R04 = N = the number of steps         R11 to R51: constants

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)
 
 
 01  LBL "T"
 02  *
 03  ST+ X
 04  CHS
 05  RTN

 
"T" ASTO 00

 0    STO 01
 1    STO 02    and if we choose  h = 0.1
0.1  STO 03
10   STO 04    XEQ "RK8"   yields            x = 1                        ( in 2mn )
                              X<>Y                     y(1) = 0.3678794412   correct to 10 D!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 28 = 256 )
 

RK8B:  RK8 for Systems of 2 Differential Equations
 
 

Data Registers:               R00:  f name         ( Registers R00 to R05 are to be initialized before executing "RK8B" )

                               R01 = x0         R04 = h = stepsize                         R06 to R17: temp    R59: counter
                               R02 = y0         R05 = N = number of steps           R18 thru  R58 = the same 41 numbers used in "RK8".
                               R03 = z0
Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)
 
 
 01  LBL "U"
 02  RCL Z
 03  *
 04  +
 05  ST+ X
 06  CHS
 07  RTN

 
"U" ASTO 00

 0    STO 01
 1    STO 02
 0    STO 03   and if we choose  h = 0.1
0.1  STO 04
10   STO 05    XEQ "RK8B"   yields        x = 1                         in 3mn10s
                              RDN                      y(1) =  0.3678794412   correct to 10D.
                              RDN                      z(1) = -0.7357588824   z-error = -10-10
 

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 28 = 256 )
 

RK8N: RK8 for Systems of n Differential Equations
 

     dy1/dx = f1(x,y1, ... ,yn)
     dy2/dx = f2(x,y1, ... ,yn)
             ..................                         with the initial values:     yi(x0)         i = 1 , ... , n

     dyn/dx = fn(x,y1, ... ,yn)
 

Data Registers:             R00 = subroutine name         ( Registers R00 to R03 & R50 to R50+n are to be initialized before executing "RK8N" )

                                        R01 = n  = number of equations = number of functions        R04 thru R13: temp        R49: temp
                                        R02 = h  = stepsize
                                        R03 = N = number of steps                                                R14 thru R45 = Constants
                                                                                                                                   R46-R47-R48 are unused
                                        R50 = x0
                                        R51 = y1(x0)                                                                           R51+n thru R50+9n: temp
                                        R52 =  y2(x0)
                                          ...............         ( initial values )

                                        R50+n = yn(x0)

Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R51+n, ... , R50+2n  respectively
                               with  x, y1, ... , yn in regiters R50, R51, ... , R50+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:

       dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1  and  N = 10
       dz/dx = x.( y + z - u )         z(0) = 1
       du/dx = x.y - z.u                u(0) = 2
 
 
 01  LBL "T"
 02  RCL 51
 03  RCL 52
 04  RCL 53
 05  *
 06  *
 07  CHS
 08  STO 54
 09  RCL 51
 10  RCL 52
 11  +
 12  RCL 53
 13  -
 14  RCL 50
 15  *
 16  STO 55
 17  LASTX
 18  RCL 51
 19  *
 20  RCL 52
 21  RCL 53
 22  *
 23  -
 24  STO 56
 25  RTN

 
 "T"  ASTO 00

  Initial values:

      0   STO 50
      1   STO 51    STO 52
      2   STO 53

 n = 3      STO 01    ( 3 functions )
 h = 0.1   STO 02
 N = 10   STO 03    XEQ "RK8N"  >>>>      x   = 1                     = R50                      ( in 9mn31s )
                                                       RDN    y(1) = 0.258207906  = R51                                                                         y(1) = 0.258207906
                                                       RDN    z(1) = 1.157623979  = R52         the exact results, rounded to 9D are          z(1) = 1.157623981
                                                       RDN    u(1) = 0.842178313  = R53                                                                         u(1) = 0.842178312

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if needed )
 

Runge-Kutta-Felhberg 4-5 method
 

-In this Runge-Kutta-Fehlberg method, a 4th-order formula is used to compute the solution
 and the difference between this 4th-order formula and a 5th-order formula provides an error estimate.
 
 

Data Registers:               R00:  f name          ( Registers R00 to R04 are to be initialized before executing "RKF" )

                               R01 = x0         R03 = h = stepsize                          R05 = error estimate
                               R02 = y0         R04 = N = number of steps             R06 thru R10: temp

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)
 
 
 01  LBL "T"
 02  *
 03  ST+ X
 04  CHS
 05  RTN

 
"T" ASTO 00

 0    STO 01
 1    STO 02    and if we choose  h = 0.1
0.1  STO 03
10   STO 04    XEQ "RKF"   yields          x = 1
                              X<>Y                    y(1) = 0.367879263

and the error estimate    RCL 05      >>>       -9.7 10-8

-The error is actually  -1.8 10-7
-The successive errors are simply added to R05, so they are sometimes underestimated...
-Press R/S to continue with the same h and N.
 

RKFB:  RKF for Systems of 2 Differential Equations
 
 

Data Registers:                R00:  f name                    ( Registers R00 to R05  are to be initialized before executing "RKFB" )

                               R01 = x0         R04 = h = stepsize                         R06 = y-error estimate           R08 to R16: temp
                               R02 = y0         R05 = N = number of steps           R07 = z-error estimate
                               R03 = z0
Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)
 
 
 01  LBL "U"
 02  RCL Z
 03  *
 04  +
 05  ST+ X
 06  CHS
 07  RTN

 
"U" ASTO 00

 0    STO 01
 1    STO 02
 0    STO 03    and if we choose  h = 0.1
0.1  STO 04
10   STO 05    XEQ "RKFB"   yields       x  =  1                      ( in 2mn16s )
                              RDN                      y(1) =  0.367879517
                              RDN                      z(1) = -0.735759034

and the error    RCL 06     >>>    -8.7 10-8          actually   y-error =    7.6 10-8
 estimates:        RCL 07     >>>    -2.1 10-7          actually    z-error =  -1.5 10-7

-Press R/S to continue with the same h and N.
 

RKFN:  RKF for Systems of n Differential Equations
 

     dy1/dx = f1(x,y1, ... ,yn)
     dy2/dx = f2(x,y1, ... ,yn)
             ..................                         with the initial values:     yi(x0)         i = 1 , ... , n

     dyn/dx = fn(x,y1, ... ,yn)
 

Data Registers:             R00 = subroutine name       ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RKFN" )

                                        R01 = n  = number of equations = number of functions         R04 thru R14: temp        R15 thru R19: unused
                                        R02 = h  = stepsize
                                        R03 = N = number of steps

                                        R20 = x0                                         at the end:
                                        R21 = y1(x0)                                     R21+n = err(y1)                                      R21+n thru R20+8n: temp
                                        R22 =  y2(x0)                                    R22+n = err(y2)
                                          ...............                                      ......................

                                        R20+n = yn(x0)                                R20+2n = err(yn)

                   Where  err(yi) is an error-estimate of yi    ( during the calculations, these error-estimates are moved in registers R21+2n .... R20+3n )

Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+n, ... , R20+2n  respectively
                           with  x, y1, ... , yn in regiters R20, R21, ... , R20+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:

       dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1  &  N = 10
       dz/dx = x.( y + z - u )         z(0) = 1
       du/dx = x.y - z.u                u(0) = 2
 
 
 01  LBL "T"
 02  RCL 21
 03  RCL 22
 04  RCL 23
 05  *
 06  *
 07  CHS
 08  STO 24
 09  RCL 21
 10  RCL 22
 11  +
 12  RCL 23
 13  -
 14  RCL 20
 15  *
 16  STO 25
 17  LASTX
 18  RCL 21
 19  *
 20  RCL 22
 21  RCL 23
 22  *
 23  -
 24  STO 26
 25  RTN

 
 "T"  ASTO 00

 Initial values:

           0   STO 20
           1   STO 21    STO 22
           2   STO 23

 n = 3      STO 01    ( 3 functions )
 h = 0.1   STO 02
 N = 10   STO 03    XEQ "RKF"  >>>>      x   = 1                     = R20                      ( in 5mn40s )
                                                     RDN    y(1) = 0.258207319  = R21                                                              err(y) = -603 E-9 = R24
                                                     RDN    z(1) = 1.157624973  = R22                and the error estimates            err(z) = 1062 E-9 = R25
                                                     RDN    u(1) = 0.842178529  = R23                                                              err(u) = 1768 E-9 = R26

-Actually, the true errors are    -587 E-9  for y(1)     992 E-9  for z(1)   217 E-9  for u(1)

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if needed )
 

RKS2:  RK4 for Systems of 2 Second order Conservative Equations
 

-"RKS2" solves the system:

     y" = f(x,y,z)         y(x0) = y0     y'(x0) = y'0
     z" = g(x,y,z)         z(x0) = z0     z'(x0) = z'0

 by the 4th order formulae:

   y(x+h) = y(x) + h ( y'(x) + k1/6 + 2k2/3 )
   y'(x+h) = y'(x) + k1/6 + 2k2/3 + k3/6                 where  k1 = h.f (x,y)  ;  k2 = h.f(x+h/2,y+h.y'/2+h.k1/8)  ;  k3 = h.f(x+h,y+h.y'+h.k2/2)
 

Data Registers:       R00:  subroutine name                         ( Registers R00 thru R07  are to be initialized before executing "RKS2" )

                                  R01 = x0         R04 = h = stepsize                            R06 = y'0          R08 to R13: temp
                                  R02 = y0         R05 = N = number of steps               R07 = z'0
                                  R03 = z0
Flags: /
Subroutine:  a program which calculates f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /     z(x0+N.h)
           Y             /     y(x0+N.h)
           X             /       x0+N.h

 
Example:

     d2y/dx2 =  -y.z                  y(0) = 2    y'(0) = 1         Evaluate  y(1) , z(1) , y'(1) , z'(1)    with  h = 0.1  &  N = 10
     d2z/dx2 =  x.( y + z )         z(0) = 1    z'(0) = 1
 
 
  LBL "T"
  RDN
  STO Z
  X<>Y
  CHS
  ST* Z
  -
  R^
  *
  RTN

 
   "T"   ASTO 00

    0     STO 01              1  STO 06
    2     STO 02              1  STO 07
    1     STO 03
   0.1   STO 04
   10    STO 05

 XEQ "RKS2"  >>>>     x   =  1                                                                                                  ( in 41 seconds )
                         RDN   y(1) = 1.531358015     and   R06 = y'(1) = -2.312838895
                         RDN   z(1) = 2.620254480              R07 = z'(1) =  2.941751649

-With  h = 0.05  it yields   y(1) =  1.531356736      y'(1) = -2.312840085
                                        z(1) =  2.620254295      z'(1) =  2.941748608

-Actually, the exact results - rounded to 9D - are

      y(1) =  1.531356646      y'(1) = -2.312840137
      z(1) =  2.620254281      z'(1) =  2.941748399

-To continue the calculations, simply press R/S  ( after changing h/2 & N in registers R04 & R05 if needed )
 

RKS3:  RK4 for Systems of 3 Second order Conservative Equations
 

-"RKS3" solves the system:

     y" = f(x,y,z,u)         y(x0) = y0     y'(x0) = y'0
     z" = g(x,y,z,u)         z(x0) = z0     z'(x0) = z'0
     u" = h(x,y,z,u)        u(x0) = u0     u'(x0) = u'0
 

Data Registers:       R00:  subroutine name                         ( Registers R00 thru R09  are to be initialized before executing "RKS3" )

                                  R01 = x0         R05 = h  = stepsize                           R07 = y'0
                                  R02 = y0         R06 = N = number of steps               R08 = z'0
                                  R03 = z0                   R10 to R17: temp                     R09 = u'0
                                  R04 = u0
Flags: /
Subroutine:  a program which calculates f(x;y;z;u) in Z-register , g(x;y;z;u) in Y-register and  h(x;y;z;u) in X-register
                       assuming x is in X-register , y is in Y-register , z is in Z-register and u is in T-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           T             /       u(x0+N.h)
           Z             /       z(x0+N.h)
           Y             /       y(x0+N.h)
           X             /         x0+N.h

 
Example:

     d2y/dx2 = -y.z.u                    y(0) = 1        y'(0) = 1         Evaluate  y(1) , z(1) , u(1)  ,   y'(1) , z'(1) , u'(1)      with  h = 0.1  &  N = 10
     d2z/dx2 = x.( y + z - u )         z(0) = 1        z'(0) = 1
     d2u/dx2 = x.y - z.u                 u(0) = 2       u'(0) = 1
 
 
  LBL "T"
  R^
  CHS
  STO L
  R^
  ST* Y
  ST+ L
  CLX
  RCL Z
  R^
  ST+ L
  ST* Y
  X<> Z
  ST+ Y
  ST* Z
  X<> T
  LASTX
  *
  X<>Y
  RTN

 
 "T"   ASTO 00

  0     STO 01
  1     STO 02  STO 03  STO 07  STO 08  STO 09
  2     STO 04
 0.1   STO 05
 10    STO 06    XEQ "RKS3"  >>>>     x   = 1                                                                                             ( in 65 seconds )
                                                 RDN   y(1) = 0.439528419           y'(1) = -2.101120400 = R07
                                                 RDN   z(1) = 2.070938499           z'(1) =   1.269599239 = R08
                                                 RDN   u(1) = 1.744522976           u'(1) = -1.704232092 = R09

-With  h = 0.05  it yields

      y(1) = 0.439524393          y'(1) = -2.101122784
      z(1) = 2.070940521          z'(1) =   1.269597110
      u(1) = 1.744524843          u'(1) = -1.704234567

-Actually, the exact results rounded to 9D are:

      y(1) = 0.439524100          y'(1) = -2.101122880
      z(1) = 2.070940654          z'(1) =   1.269596950
      u(1) = 1.744524964          u'(1) = -1.704234756

-Press R/S  to continue the calculations ( after changing h & N in registers R05 & R06 if needed )
-The subroutine may be difficult to write, but "RKS3" is faster than "RKSN" below, for 3 equations.
 

RKSN:  RK4 for Systems of n Second order Conservative Equations
 

-"RKSN"  solves the system

  y"1 = f1 ( x, y1 , ... , yn)               y1(x0) = y1,0     y'1(x0) = y'1,0
  ..................................                 ........................................

  y"n = fn ( x, y1 , ... , yn)               yn(x0) = yn,0     y'n(x0) = y'n,0
 

Data Registers:             R00 = subroutine name       ( Registers R00 thru R03 and R10 thru R10+2n are to be initialized before executing "RKSN" )

                                        R01 = n  = number of equations = number of functions         R04 thru R09: temp
                                        R02 = h  = stepsize
                                        R03 = N = number of steps

                                        R10 = x0
                                        R11 = y1(x0)                                    R11+n = y'1(x0)                                      R11+2n thru R10+6n: temp
                                        R12 =  y2(x0)                                   R12+n = y'2(x0)
                                          ...............                                      ......................

                                        R10+n = yn(x0)                                R10+2n = y'n(x0)

                                      ( during the calculations, y'i  are moved in registers R11+2n .... R10+3n )

Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R11+n, ... , R10+2n  respectively
                           with  x, y1, ... , yn in regiters R10, R11, ... , R10+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:

    d2y/dx2 = -y.z.u                    y(0) = 1        y'(0) = 1         Evaluate  y(1) , z(1) , u(1)  ,   y'(1) , z'(1) , u'(1)      with  h = 0.1  &  N = 10
    d2z/dx2 = x.( y + z - u )         z(0) = 1        z'(0) = 1
    d2u/dx2 = x.y - z.u                 u(0) = 2       u'(0) = 1
 
 
 01  LBL "T"
 02  RCL 11
 03  RCL 12
 04  RCL 13
 05  *
 06  *
 07  CHS
 08  STO 14
 09  RCL 11
 10  RCL 12
 11  +
 12  RCL 13
 13  -
 14  RCL 10
 15  *
 16  STO 15
 17  RCL 10
 18  RCL 11
 19  *
 20  RCL 12
 21  RCL 13
 22  *
 23  -
 24  STO 16
 25  RTN

 
 "T"  ASTO 00

 Initial values:     0   STO 10

                         1   STO 11    STO 12    STO 14    STO 15    STO 16
                         2   STO 13

 n = 3      STO 01    ( 3 functions )
 h = 0.1   STO 02
 N = 10   STO 03    XEQ "RKSN"  >>>>    x   =  1                     = R10                                                                ( in 2mn19s )
                                                       RDN   y(1) = 0.439528419  = R11          y'(1) = -2.101120401 = R14
                                                       RDN   z(1) = 2.070938499  = R12          z'(1) =   1.269599239 = R15
                                                       RDN   u(1) = 1.744522977  = R13          u'(1) = -1.704232091 = R16
 

-Press R/S  to continue the calculations ( after changing h & N in registers R02 & R03 if need be )
 

Runge-Kutta-Verner 5-6 method
 
 

-All n stage explicit Runge-Kutta formulae with error estimate are determined by ( n2 + 5.n - 2 ) / 2 coefficients ai,j ; bi ; ci ; di

    k1 = h.f(x;y)
    k2 = h.f(x+b2.h;y+a2,1.k1)
    k3 = h.f(x+b3.h;y+a3,1.k1+a3,2.k2)
    .............................................................                                         (  actually,     ai,1 + ai,2 + .......... + ai,i-1 = bi  )

    kn = h.f(x+bn.h;y+an,1.k1+an,2.k2+ .......... +an,n-1.kn-1)

then,     y(x+h) = y(x) + c1.k1+ ................ + cn.kn    and the difference between the higher-order and the lower-order formulae gives an error estimate:

         err-estim =  d1.k1+ ................ + dn.kn

-"RKV" uses an 8-stage method with order 5 with a 6th-order error-estimating formula.
 

Data Registers:            R00:  f name                      ( Registers R00 to R05 are to be initialized before executing "RKV" )

                               R01 = x0                                    R05 = y-error                       R19 to R69 =  the 51 coefficients given by "C51"
                               R02 = y0                                      R06 to R10: temp
                               R03 = h = stepsize                        R11 to R18 = ki
                               R04 = N = number of steps

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:   Our usual equation:   y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)  with  h = 0.1  and  N = 10
 
 
 01  LBL "T"
 02  *
 03  ST+ X
 04  CHS
 05  RTN

 
"T" ASTO 00

 0    STO 01
 1    STO 02
0.1  STO 03
10   STO 04    XEQ "RKV"   yields          x = 1                                    ( in 3mn56s )
                              X<>Y                     y(1) = 0.367879457

  and the error estimate         >>>          R05 = -13 10-9

-The error is actually  16 10-9
-Press R/S to continue ( after changing  h and N in registers R03 and R04 if need be ).
 

RKVB:  RKV for Systems of 2 Differential Equations
 

Data Registers:           R00:  f name                       ( Registers R00 to R06  are to be initialized before executing "RKVB" )

                               R01 = x0                                       R06 = y-error               R30 to R80 =  the 51 coefficients given by "C51"
                               R02 = y0                                       R07 = z-error
                               R03 = y0
                               R04 = h =  stepsize                        R08  to  R13: temp
                               R05 = N = number of steps           R14    to R21 = ki(y)       R22 to R29 = ki(z)

Flags: /
Subroutine:    a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                        assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)  with  h = 0.1  and  N = 10
 
 
 01  LBL "U"
 02  RCL Z
 03  *
 04  +
 05  ST+ X
 06  CHS
 07  RTN

 
"U" ASTO 00

 0    STO 01
 1    STO 02
 0    STO 03
0.1  STO 04
10   STO 05   XEQ "RKVB"   yields       x  =  1                      ( in 6mn34s )
                              RDN                      y(1) =  0.367879378
                              RDN                      z(1) = -0.735758757

and the error estimates    RCL 06  =    -85 10-9       actually    y-error =   -63 10-9
                                      RCL 07  =    153 10-9       actually    z-error =   126 10-9

-Press R/S to continue ( after changing  h and N in registers R04 and R05 if needed ).
 

Bulirsh-Stoer Method
 

-"BST" solves the differential equation   dy/dx = f(x,y)   with the initial value  y(x0) = y0
-It uses an extrapolation to the limit and the modified midpoint formula to compute y(x0+H)
-Let

                       z0 = y0
                       z1 = z0 + h.f(x0,z0)
                       z2 = z0 + 2h.f(x0+h,z1)
                       z3 = z1 + 2h.f(x0+2h,z2)                                     where  h = H/n
                       ..................................
                       zn = zn-2 + 2h.f(x0+(n-1).h,zn-1)

    y(x0+H) ~  yn = (1/2).[ zn + zn-1 + h.f(x0+H,zn)

-The numbers yn are computed for n = 2 , 4 , 6 , ..... , 16  ( at most )  and the Lagrange extrapolation polynomial is used
  until the estimated error is smaller than a specified tolerance ( tol )

-If the estimated error is still greater than tol with  n = 16 ,  H is halved
-On the other hand, if the desired accuracy is satisfied with n < 8 , H is doubled
-Other ( and much better ) methods for stepsize control may be used ( cf "Numerical Recipes" )

-The Bulirsch-Stoer technique is quite similar to the Romberg integration:
  the modified midpoint formula is second-order only,
  but, for instance, with n = 16 and the deferred approach to the limit, we actually use a 16th-order method !
 

Data Registers:             R00 = f name                                                 ( Registers R00 thru R03 are to be initialized before executing "BST" )

                                        R01 = x0      R04 = x1     R07 = h = H/n             R13 = next to last H
                                        R02 = y0      R05 = H     R08 to R12: temp        R14, R15, ........ , R21 = y-approximations
                                        R03 = tol     R06 = n = 2 , 4 , .... , 16

                          >>>>   where tol is a small positive number that defines the desired accuracy

Flags:  F09 - F10
Subroutine:  A program that computes  f(x,y)  assuming x is in X-register and y is in Y-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Y             /          y(x1)
           X            x1            x1

 
Example:     dy/dx =   x.(y/2)2           y(0) = 1         Evaluate  y(2)  and  y(2.5)
 
 
  LBL "T"
  X<>Y 
  2
  /
  X^2
  *
  RTN

 
 "T"  ASTO 00

  0  STO 01
  1  STO 02   and if we choose  tol = 10 -7   STO 03

  2  XEQ "BST"  >>>>        x   =  2                              ( in 1mn49s )                             y(1)  has been computed with  H = 1 , n = 10
                           RDN      y(2) =  2.000000018      the exact result is 2                         y(2)  ------------------------  H = 1 , n = 14

  2.5  R/S    >>>>       x    = 2.5                                   ( in 3mn17s )
                   RDN   y(2.5) = 4.571428682     the last 3 decimals should be  571            H = 1 has been replaced by  0.5  but the tolerance 10 -7  cannot be
                                                                      the error is 1.11 10 -7                               satisfied with this stepsize. So, H finally became 0.25 ( and n = 14 )

-In fact, the solution is  y = 1/(1-x2/8)   so we'd need to increase tol in R03 and smaller and smaller stepsizes as x get closer to sqrt(8)

Note:

-Do not choose a too small tolerance , it could produce  an infinite loop.
 

BST2:  BST  for Systems of 2 Differential Equations
 

-The same method is employed to solve the system

     dy/dx = f(x,y,z)
     dz/dx = g(x,y,z)       with initial values    y(x0) = y0  ,    z(x0) = z0
 

Data Registers:             R00 = f name                                       ( Registers R00 thru R04 are to be initialized before executing "BST2" )

                                        R01 = x0      R05 = H                              R09 to R17: temp
                                        R02 = y0      R06 = x1                              R18 = next to last H-value
                                        R03 = z0      R07 = n = 2 , 4 , .... , 16      R19, R20, ........ , R26 = y-approximations
                                        R04 = tol     R08 = h = H/n                      R27, R28, ........ , R34 = z-approximations

                             >>>>   where tol is a small positive number that defines the desired accuracy

Flags:  F09 - F10
Subroutine:  A program that computes  f(x,y,z) in Y-register
                                                       and  g(x,y,z) in Z-register   assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /          z(x1)
           Y             /          y(x1)
           X            x1            x1

 
Example:    d2y/dx2 = -2y - 2x.dy/dx    y(0) = 1    y'(0) = 0       Calculate  y(1) and y'(1)        [ the exact solution is  y = exp ( -x2 ) ]

  This second-order equation is equivalent to the system

    dy/dx = z                      y(0) = 1
    dz/dx = -2y - 2x.z         z(0) = 0
 
 
  LBL "T" 
  RCL Z 
  * 
  +
  ST+ X
  CHS
  RTN

 
 "T"  ASTO 00

  0  STO 01
  1  STO 02
  0  STO 03    and with  tol = 10 -7   STO 04

  1  XEQ "BST2"  >>>>        x   =  1                                     ( in 2mn08s )
                RDN                   y(1) =  0.367879446         the last decimal should be a 1
                RDN       z(1) =  y'(1) = -0.735758909        the last 3 decimals should be 882      z-error ~ 27 E-9  <  E-7
 

BSTN:  BST  for Systems of n Differential Equations
 

     dy1/dx = f1(x,y1, ... ,yn)
     dy2/dx = f2(x,y1, ... ,yn)
             ..................                         with the initial values:     yi(x0)         i = 1 , ... , n

     dyn/dx = fn(x,y1, ... ,yn)
 

Data Registers:             R00 = subroutine name       ( Registers R00-R01-R02 and R20 thru R20+n are to be initialized before executing "BSTN" )

                                        R01 = n  = number of equations = number of functions         R03 thru R18: temp        R19: unused

                                        R02 = tolerance = a small positive number that specifies the desired accuracy

                                        R20 = x0
                                        R21 = y1(x0)                                                                           R21+n thru R20+14n: temp
                                        R22 =  y2(x0)
                                          ...............         ( initial values )

                                        R20+n = yn(x0)
Flags:  F09-F10
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+n, ... , R20+2n  respectively
                           with  x, y1, ... , yn in regiters R20, R21, ... , R20+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /         y3(x1)
           Z             /         y2(x1)
           Y             /         y1(x1)
           X            x1            x1

 
Example:

       dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  a tolerance = 10 -7
       dz/dx = x.( y + z - u )         z(0) = 1
       du/dx = x.y - z.u                u(0) = 2
 
 
 01  LBL "T"
 02  RCL 21
 03  RCL 22
 04  RCL 23
 05  *
 06  *
 07  CHS
 08  STO 24
 09  RCL 21
 10  RCL 22
 11  +
 12  RCL 23
 13  -
 14  RCL 20
 15  *
 16  STO 25
 17  LASTX
 18  RCL 21
 19  *
 20  RCL 22
 21  RCL 23
 22  *
 23  -
 24  STO 26
 25  RTN

 
 "T"  ASTO 00

 Initial values:     0   STO 20

                         1   STO 21    STO 22
                         2   STO 23

  n  =  3     STO 01    ( 3 functions )
 tol = E-7  STO 02

    1   XEQ "BSTN"  >>>>    x   =  1                    = R20                      ( in 12mn06s )
                              RDN    y(1) = 0.258207909  = R21                                                                         y(1) = 0.258207906
                              RDN    z(1) = 1.157623986  = R22         the exact results, rounded to 9D are          z(1) = 1.157623981
                              RDN    u(1) = 0.842178304  = R23                                                                         u(1) = 0.842178312
 

-The long execution time is due to the fact that the method is first used with H = 1
  but the desired accuracy is not satisfied. So the stepsize is divided by 2.

-To compute  y(2)  z(2)  u(2)  key in  2  R/S   it yields              ( in 6m46s )

                y(2) = 0.106363294                                                 y(2) = 0.106363329
                z(2) = 3.886706181         the exact results are          z(2) = 3.886706159         ( rounded to 9D )
                u(2) = 0.196515847                                                 u(2) = 0.196515847

Notes:

-The initial stepsize is always +/-1 ( line 03 )
-The modified midpoint rule ( and Richarson extrapolation ) are used with  h/2 , h/4 , h/6 , ..... , h/16 ( at most )  until the desired accuracy is satisfied.
 

Stoermer's rule for Conservative 2nd-order Equations
 

-"STRM" solves the differential equation     d2y/dx2 = f(x,y)     with initial values    y(x0) = y0  ,   y'(x0) = y'0
  where the first derivative does not appear explicitly.

-This kind of equations may be re-written as a system which could be solved by "BST2" but the following method is faster.

-Like "BST", "STRM" also uses the deferred approach to the limit but - instead of the modified midpoint method - the Stoermer's rule is employed:
-Let

    d0 = h.[ y'0 + (h/2).f(x0,y0) ]
    y1 = y0 + d0
    dk = dk-1 + h2f(x0+k.h,yk)                        dk = yk+1 - yk          k = 1 , 2 , ..... , n-1         h = H/n
    y'n = dn-1/h + (h/2).f(x0+H,yn)
 

Data Registers:             R00 = f name                                       ( Registers R00 thru R04 are to be initialized before executing "STRM" )

                                        R01 = x0      R05 = H                              R09 to R14: temp
                                        R02 = y0      R06 = x1                              R15 = next to last H-value
                                        R03 = y'0     R07 = n = 2 , 4 , .... , 16      R16, R17, ........ , R23 = y-approximations
                                        R04 = tol     R08 = h = H/n                      R24, R25, ........ , R31 = y'-approximations

                          >>>>   where tol is a small positive number that defines the desired accuracy

Flags:  F09 - F10
Subroutine:  A program that computes  f(x,y)   assuming x is in X-register and y is in Y-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /          y'(x1)
           Y             /          y(x1)
           X            x1            x1

 
Example:        d2y/dx2 = -y.( x2 + y2 )1/2         y(0) = 1  ,  y'(0) = 0       Evaluate  y(1) and y(PI)
 
 
  LBL "T"
  X^2 
  RCL Y
  X^2
  +
  SQRT
  *
  CHS
  RTN

 
"T"  ASTO 00

  0  STO 01
  1  STO 02
  0  STO 03    and with  tol = 10 -7   STO 04

  1  XEQ "STRM"  >>>>        x   =  1                                  ( in 53 seconds )
                              RDN      y(1) =  0.536630616         the 9 decimals are correct
                              RDN     y'(1) = -0.860171925         the last decimal should be a 7

 PI        R/S         >>>>        x   =  3.141592654                  ( in 2mn24s )
                            RDN    y(PI) = -0.411893053         the 9 decimals are exact
                            RDN    y'(PI) =  1.018399901         the last decimal should be a 3
 

2nd order differential equations ( 4th order method )
 

-"SRK4" solves the equation  y" = f(x,y,y')     with the initial values:   y(x0) = y0  ,   y'(x0) = y'0
 

Data Registers:             R00 = function name                        ( Registers R00 thru R05 are to be initialized before executing "SRK4"  )

                                        R01 = x0
                                        R02 = y0          R04 =  h = stepsize
                                        R03 = y'0         R05 =  m  = number of steps             R06 thru R10: temp
Flags: /
Subroutine:       A program which computes  y" = f(x,y,y')  assuming  x , y , y'  are in registers X , Y , Z ( respectively )  upon entry
 
 
      STACK        INPUTS      OUTPUTS
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:     Let's consider the Lane-Emden equation of index 3     y" = -(2/x) y' - y3   with the initial values  y(0) = 1 , y'(0) = 0

-There is already a difficulty because x = 0 is a singular point!
-But if we use a series expansion, we find after substitution that  y = 1 + a.x2 + ....  will satisfy the LEE  if   a = -1/6    whence  y"(0) = -1/3

-So, we can load the following subroutine into program memory:
 
 
  LBL "T"
  X=0?
  GTO 00
  RCL Z
  ST+ X
  X<>Y
  /
  X<>Y
  3
  Y^X
  +
  CHS
  RTN
  LBL 00
  3
  1/X
  CHS
  RTN

 
-If we want to estimate y(1)  with a stepsize h = 0.1   ( whence  m = 10 )

  "T"  ASTO 00

   0   STO 01  STO 03             0.1   STO 04  = h
   1   STO 02                            10    STO 05

  XEQ "SRK4"   >>>>       x   =  1                            ( in 56 seconds )
                          RDN     y(1)  =  0.855057170
                          RDN     y'(1) = -0.252129561

-These new "initial" values are also stored in registers R01 R02 R03
-With  h = 0.05  and  m = 20 we would have found  y(1) = 0.855057539  &  y'(1) = -0.252129276
 

Notes:    The solution of the Lane-Emden Equation of index 3 looks like this:

      y
    1|*                I
      |                  *
      |                              *
      |                                           *
      |                                                              *
      |-----------------------------------------------------------*x1-------- x
     0
              y(x1) = 0  for  x1 = 6.896848619    and    y'(x1) = -0.0424297576
and there is an inflexion point I with  xI = 1.495999168 , yI = 0.720621687 and  y'(xI) = -0.279913175

-The solutions of the Lane-Emden Equations of index n   ( y" + (2/x).y' + yn = 0  ,  y(0)= 1 , y'(0) = 0 )
  can be expressed by elementary functions for only 3 n-values:

   a)  n = 0   y(x) = 1 - x2/6
   b)  n = 1   y(x) = (sin x)/x
   c)  n = 5   y(x) = ( 1+x2/3) -1/2
 

3rd order differential equations ( 4th order method )
 

-Likewise, "TRK4" solves the equation  y"' = f(x,y,y',y")     with the initial values:   y(x0) = y0  ,   y'(x0) = y'0  ,  y"(x0) = y"0

Data Registers:             R00 = function name                          ( Registers R00 thru R06 are to be initialized before executing "TRK4"  )

                                        R01 = x0
                                        R02 = y0
                                        R03 = y'0         R05 =  h  = stepsize
                                        R04 = y"0        R06 =  m  = number of steps             R07 thru R13: temp
Flags: /
Subroutine:       A program which computes  y"' = f(x,y,y',y")  assuming  x , y , y' , y"  are in registers X , Y , Z , T ( respectively )  upon entry
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y"(x0+m.h)
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:      y"' =  2xy" - x2y' + y2   with  y(0) = 1 , y'(0) = 0 , y"(0) = -1

-We load, for instance:
 
 
  LBL "S"
  X^2
  ST* Z
  X<> L
  ST+ X
  ST* T
  RDN
  X^2
  -
  -
  RTN

 
-With  h =  0.1  and  m = 10

  "S"  ASTO 00

   0   STO 01  STO 03                  0.1    STO 05 = h
   1   STO 02  CHS  STO 04         10    STO 06

  XEQ "TRK4"   >>>>       x   =  1                            ( in 49 seconds )
                          RDN     y(1)  =  0.595434736
                          RDN     y'(1) = -0.776441445
                          RDN    y"(1) =  -0.791715205

-With  h = 0.05  and  m = 20,  we would find  y(1) = 0.595431304  ,  y'(1) = -0.776444326  ,  y"(1) = -0.791718300
 

N-th order differential equations ( 4th order method )
 

-The differential equation is now  y(n) = f(x,y,y',y",.....,y(n-1))     with the initial values:   y(x0) = y0  ,   y'(x0) = y'0  ,  ........  ,  y(n-1)(x0) = y(n-1)0
 

Data Registers:             R00 = function name           ( Registers R00 thru R03 and R09 thru Rn+9 are to be initialized before executing "NRK4"  )

                                        R01 =  n      ( order )
                                        R02 =  h =  stepsize                                        R04 thru R08 & Rn+10 thru R3n+9: temp
                                        R03 =  m  =  number of steps

                                        R09 = x0       R10 = y0       R11 = y'0      R12 = y"0   .....................    Rn+9 = y(n-1)0
Flags: /
Subroutine:   A program which computes   y(n) = f(x,y,y',y",.....,y(n-1))   assuming  x , y , y' , y" , ........ , y(n-1)  are in registers  R09 R10 R11 R12 .... Rnn+9
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y"(x0+m.h)
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:      y(5) = y(4) - 2x.y"' + y" - y.y'   with  y(0) = 1 , y'(0) = y"'(0) = y(4)(0) = 0 , y"(0) = -1

-We load, for instance:
 
 
  LBL "W"
  RCL 14
  RCL 09
  ST+ X
  RCL 13
  *
  -
  RCL 12
  +
  RCL 10
  RCL 11
  *
  -
  RTN

 
-With  h =  0.1  and  m = 10

  "W"   ASTO 00                                                       and the initial values:     0  STO 09   STO 11  STO 13  STO 14
                                                                                                                     1  STO 10  CHS  STO 12
     5    STO 01        ( fifth order equation )
   0.1   STO 02        ( h )
   10    STO 03        ( m = 10 )

  XEQ "NRK4"   >>>>       x   =  1                      = R09                  ( in 2mn48s )
                          RDN     y(1)  =  0.491724880   = R10
                          RDN     y'(1) = -1.041200697   = R11       and     RCL 13 =  y"'(1)  = -0.479803795
                          RDN    y"(1) =  -1.163353624  = R12                  RCL 14 = y(4)(1) = -0.897595630

-With  h = 0.05  and  m = 20,  it yields:

     y(1) = 0.491724223 ,  y'(1) = -1.041200398 ,  y"(1) = -1.163353549 ,  y"'(1)  = -0.479804004 ,  y(4)(1) = -0.897594479
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4
[3]  F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[4]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5