hp41programs

Diffeq Systems of Differential Equations for the HP-41
 

Overview
 

 1°) Systems of first-order Differential Equations

   a) The "Classical" 4th-order Runge-Kutta Method
   b) A Runge-Kutta Method of order 6
   c) A Runge-Kutta Method of order 8
   d) Runge-Kutta-Fehlberg 4-5 Method
   e) Bulirsch-Stoer Method

 2°) Systems of second-order Conservative Equations

  a)  Systems of 2 Equations
  b)  Systems of 3 Equations
  c)  Systems of n Equations
 

-The programs listed in paragraph 1°)  solve a system of n first-order 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)

-In paragraph 2°) the HP-41 solves systems of second-order equations where the first derivatives  y'i(x) do not appear explicitly:

     d2y1/dx2 = f1(x,y1, ... ,yn)
     d2y2/dx2 = f2(x,y1, ... ,yn)
             ..................                                        with the initial values:     yi(x0)   ,    y'i(x0)       i = 1 , ... , n

     d2yn/dx2 = fn(x,y1, ... ,yn)

-All these programs use 2 instructions from the X-functions module:  REGMOVE  and  REGSWAP - Except in §2°)a)b)
 

1°) Systems of first-order Differential Equations
 

     a) The "Classical" 4th-order Runge-Kutta Method
 
 

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

                                      •  R01 = n  = number of equations = number of functions         R04 thru R09: temp
                                      •  R02 = h/2  = stepsize / 2
                                      •  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
 

-If you don't have an HP-41 CX, line 24 can be replaced by the 5 lines:

    0
    LBL 00
    STO IND Y
    ISG Y
    GTO 00
 
 

01  LBL "RK4" 
02  RCL 03
03  STO 09
04  RCL 01
05  RCL 01
06  RCL 01
07  10.01
08  +
09  STO 04
10  INT
11  +
12  STO 05
13  +
14  STO 06
15  +
16  STO 07
17  RCL 05
18  RCL 06
19  1
20  +		 21   E3
22  /
23  +
24  CLRGX
25  11
26  LASTX
27  +
28  RCL 01
29   E6
30  /
31  +
32  STO 08
33  GTO 03
34  LBL 01
35  XEQ IND 00
36  LBL 02
37  RCL IND 05
38  RCL 02
39  *
40  ST+ IND 06		 41  FS? 06
42  ST+ IND 06
43  FC? 07
44  ST+ X
45  RCL IND 07
46  +
47  STO IND 04
48  DSE 07
49  DSE 06
50  DSE 05
51  DSE 04
52  GTO 02
53  RCL 01
54  ST+ 04
55  ST+ 05
56  ST+ 06
57  ST+ 07
58  RTN
59  LBL 03
60  RCL 08		 61  REGMOVE
62  CF 06
63  SF 07
64  XEQ 01
65  SF 06
66  RCL 02
67  ST+ 10
68  XEQ 01
69  CF 07
70  XEQ 01
71  CF 06
72  RCL 02
73  ST+ 10
74  XEQ 01
75  RCL 08
76  REGSWAP
77  RCL 04
78  RCL 06
79  3
80  SIGN		 81  LBL 04
82  CLX
83  X<> IND Y
84  LASTX
85  /
86  ST+ IND Z
87  DSE Y
88  DSE Z
89  GTO 04
90  DSE 09
91  GTO 03
92  RCL 13
93  RCL 12
94  RCL 11
95  RCL 10
96  END
 

 
   ( 155 bytes / SIZE 4n+11 )
 
 

      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  STO16
 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/2 is stored in R02                0.05  STO 02
   N is stored in R03                 10    STO 03   and  XEQ "RK4"

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,
but naturally, roundoff errors ( and execution time ) will limit the accuracy attainable!
 
 

     b) A Runge-Kutta Method of Order 6
 
 

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

                                      •  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

-Line 301 is a three-byte  GTO 01
 
 

  01  LBL "RK6"
  02  RCL 03
  03  STO 04
  04  21.021
  05  RCL 01
  06   E6
  07  /
  08  +
  09  RCL 01
  10  7
  11  *
  12  +
  13  STO 05
  14  RCL 01
  15  RCL 01
  16  RCL 01
  17  20.02
  18  +
  19  STO 06
  20  +
  21  STO 07
  22  +
  23  STO 08
  24  +
  25  STO 09
  26  +
  27  STO 10
  28  +
  29  STO 11
  30  +
  31  STO 12
  32  LBL 01
  33  RCL 20
  34  STO 13
  35  RCL 05
  36  REGSWAP
  37  REGMOVE
  38  XEQ IND 00
  39  LBL 02
  40  RCL IND 07
  41  RCL 02
  42  *
  43  STO IND 08
  44  3
  45  /
  46  ST+ IND 06
  47  DSE 08
  48  DSE 07
  49  DSE 06
  50  GTO 02
  51  RCL 01		   52  ST+ 06
  53  ST+ 07
  54  ST+ 08
  55  RCL 02
  56  LASTX
  57  /
  58  ST+ 20
  59  XEQ IND 00
  60  RCL 05
  61  REGMOVE
  62  LBL 03
  63  RCL IND 07
  64  RCL 02
  65  *
  66  STO IND 09
  67  1.5
  68  /
  69  ST+ IND 06
  70  DSE 09
  71  DSE 07
  72  DSE 06
  73  GTO 03
  74  RCL 01
  75  ST+ 06
  76  ST+ 07
  77  ST+ 09
  78  RCL 02
  79  3
  80  /
  81  ST+ 20
  82  XEQ IND 00
  83  RCL 05
  84  REGMOVE
  85  LBL 04
  86  RCL IND 07
  87  RCL 02
  88  *
  89  STO IND 10
  90  RCL IND 09
  91  4
  92  *
  93  -
  94  RCL IND 08
  95  -
  96  12
  97  /
  98  ST- IND 06
  99  DSE 10
100  DSE 09
101  DSE 08
102  DSE 07		 103  DSE 06
104  GTO 04
105  RCL 01
106  ST+ 06
107  ST+ 07
108  ST+ 08
109  ST+ 09
110  ST+ 10
111  RCL 02
112  3
113  /
114  ST- 20
115  XEQ IND 00
116  RCL 05
117  REGMOVE
118  LBL 05
119  RCL IND 07
120  RCL 02
121  *
122  STO IND 11
123  18
124  *
125  RCL IND 10
126  7
127  *
128  +
129  RCL IND 09
130  22
131  *
132  -
133  RCL IND 08
134  5
135  *
136  +
137  9.6
138  /
139  ST+ IND 06
140  DSE 11
141  DSE 10
142  DSE 09
143  DSE 08
144  DSE 07
145  DSE 06
146  GTO 05
147  RCL 01
148  ST+ 06
149  ST+ 07
150  ST+ 08
151  ST+ 09
152  ST+ 10
153  ST+ 11		 154  RCL 02
155  2
156  /
157  ST+ 20
158  XEQ IND 00
159  RCL 05
160  REGMOVE
161  LBL 06
162  RCL IND 07
163  RCL 02
164  *
165  STO IND 12
166  5
167  /
168  RCL IND 11
169  +
170  RCL IND 10
171  4
172  /
173  -
174  RCL IND 08
175  18
176  *
177  RCL IND 09
178  55
179  *
180  -
181  60
182  /
183  +
184  2
185  /
186  ST+ IND 06
187  DSE 12
188  DSE 11
189  DSE 10
190  DSE 09
191  DSE 08
192  DSE 07
193  DSE 06
194  GTO 06
195  RCL 01
196  ST+ 06
197  ST+ 07
198  ST+ 08
199  ST+ 09
200  ST+ 10
201  ST+ 11
202  ST+ 12
203  RCL 02
204  6		 205  /
206  RCL 13
207  +
208  STO 20
209  XEQ IND 00
210  RCL 05
211  REGMOVE
212  LBL 07
213  RCL IND 07
214  RCL 02
215  *
216  80
217  X<>Y
218  ST* Y
219  X<> IND 09
220  99
221  *
222  +
223  RCL IND 11
224  118
225  *
226  -
227  20
228  *
229  RCL IND 12
230  ST+ IND 09
231  128
232  *
233  +
234  RCL IND 10
235  ST+ IND 11
236  215
237  *
238  +
239  RCL IND 08
240  783
241  *
242  -
243  780
244  /
245  ST+ IND 06
246  DSE 12
247  DSE 11
248  DSE 10
249  DSE 09
250  DSE 08
251  DSE 07
252  DSE 06
253  GTO 07
254  RCL 01
255  ST+ 06		 256  ST+ 07
257  ST+ 08
258  ST+ 09
259  ST+ 10
260  ST+ 11
261  ST+ 12
262  RCL 02
263  RCL 13
264  +
265  STO 20
266  XEQ IND 00
267  RCL 05
268  REGMOVE
269  LBL 08
270  RCL IND 07
271  RCL 02
272  *
273  RCL IND 08
274  +
275  13
276  *
277  RCL IND 09
278  32
279  *
280  +
281  RCL IND 11
282  55
283  *
284  +
285  .5
286  %
287  ST+ IND 06
288  DSE 11
289  DSE 09
290  DSE 08
291  DSE 07
292  DSE 06
293  GTO 08
294  RCL 01
295  ST+ 06
296  ST+ 07
297  ST+ 08
298  ST+ 09
299  ST+ 11
300  DSE 04
301  GTO 01
302  RCL 23
303  RCL 22
304  RCL 21
305  RCL 20
306  END

 
   ( 498 bytes / SIZE 21+8n )
 
 

      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
 

  LBL "T"      *                 RCL 21      -                   RCL 20      RCL 23       RTN            "T"  ASTO 00
  RCL 21      *                 RCL 22      RCL 20        RCL 21      *
  RCL 22      CHS           +                 *                  *                 -
  RCL 23      STO 24      RCL 23       STO 25        RCL 22      STO 26

 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 "RK6"  >>>>      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 )
 
 

     c) A Runge-Kutta Method of Order 8
 
 

Data Registers:           •  R00 = subroutine name     ( Registers R00 to R03 , R14 to R48 and R50 to R50+n are to be initialized before executing "RK8" )

                                      •  R01 = n  = number of equations = number of functions        R04 thru R13: temp        R49: temp
                                      •  R02 = h  = stepsize
                                      •  R03 = N = number of steps                                            •  R14 thru R48 = Constants

                                      •  R50 = x0
                                      •  R51 = y1(x0)                                                                           R51+n thru R50+9n: temp
                                      •  R52 =  y2(x0)
                                          ...............         ( initial values )

                                      •  R50+n = yn(x0)
 

           •  R14 =  0.8961811934                 •  R23 =  0.1996369936                    •  R32 =  0.9576053440                    •  R41 =  0.9401094520
           •  R15 = -0.2117115009                •  R24 =  0.09790045616                  •  R33 =  -0.6325461607                   •  R42 = -2.229158210
           •  R16 =  0.8273268354                 •  R25 =  0.3285172131                    •  R34 =  0.1527777778                    •  R43 =  7.553840442
           •  R17 =  0.06514850915               •  R26 = -0.3497052863                    •  R35 = -0.009086961101               •  R44 = -7.164951553
           •  R18 =  0.5766714727                 •  R27 = -0.03302551131                  •  R36 =  1.062498447                      •  R45 =  2.451380432
           •  R19 =  0.1855068535                 •  R28 =  0.1289862930                     •  R37 = -1.810863083                     •  R46 =  0.05
           •  R20 =  0.3865246267                 •  R29 =  0.1726731646                     •  R38 =  2.031083139                     •  R47 =   0.2722222222
           •  R21 = -0.4634553896                •  R30 =  -0.01186868389                  •  R39 = -0.6379313502                   •  R48 =  0.3555555556
           •  R22 =  0.3772937693                 •  R31 =  0.002002165993                 •  R40 =  -0.5512205631
 

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

-Line 459 is a three-byte GTO 01
 
 

  01  LBL "RK8"
  02  51.051
  03  RCL 01
  04   E6
  05  /
  06  +
  07  RCL 01
  08  8
  09  *
  10  +
  11  STO 05
  12  RCL 01
  13  RCL 01
  14  RCL 01
  15  50.05
  16  +
  17  STO 06
  18  +
  19  STO 07
  20  +
  21  STO 08
  22  +
  23  STO 09
  24  +
  25  STO 10
  26  +
  27  STO 11
  28  +
  29  STO 12
  30  +
  31  STO 13
  32  RCL 03
  33  STO 49
  34  GTO 01
  35  LBL 00
  36  XEQ IND 00
  37  RCL 05
  38  REGMOVE
  39  RTN
  40  LBL 01
  41  RCL 50
  42  STO 04
  43  RCL 05
  44  REGSWAP
  45  REGMOVE
  46  XEQ IND 00
  47  LBL 02
  48  RCL IND 07
  49  RCL 02
  50  *
  51  STO IND 08
  52  2
  53  /
  54  ST+ IND 06
  55  DSE 08
  56  DSE 07
  57  DSE 06
  58  GTO 02
  59  RCL 01
  60  ST+ 06
  61  ST+ 07
  62  ST+ 08
  63  RCL 02
  64  LASTX
  65  /
  66  ST+ 50
  67  XEQ 00
  68  LBL 03
  69  RCL IND 07
  70  RCL 02
  71  *
  72  STO IND 09
  73  RCL IND 08
  74  +
  75  4
  76  /
  77  ST+ IND 06
  78  DSE 09		   79  DSE 08
  80  DSE 07
  81  DSE 06
  82  GTO 03
  83  RCL 01
  84  ST+ 06
  85  ST+ 07
  86  ST+ 08
  87  ST+ 09
  88  XEQ 00
  89  LBL 04
  90  RCL IND 07
  91  RCL 02
  92  *
  93  STO IND 10
  94  RCL 14
  95  *
  96  RCL IND 09
  97  RCL 15
  98  *
  99  +
100  RCL IND 08
101  7
102  /
103  +
104  ST+ IND 06
105  DSE 10
106  DSE 09
107  DSE 08
108  DSE 07
109  DSE 06
110  GTO 04
111  RCL 01
112  ST+ 06
113  ST+ 07
114  ST+ 08
115  ST+ 09
116  ST+ 10
117  RCL 02
118  RCL 16
119  *
120  RCL 04
121  +
122  STO 50
123  XEQ 00
124  LBL 05
125  RCL IND 07
126  RCL 02
127  *
128  STO IND 11
129  RCL 17
130  *
131  RCL IND 10
132  RCL 18
133  *
134  +
135  RCL IND 08
136  RCL 19
137  *
138  +
139  ST+ IND 06
140  DSE 11
141  DSE 10
142  DSE 08
143  DSE 07
144  DSE 06
145  GTO 05
146  RCL 01
147  ST+ 06
148  ST+ 07
149  ST+ 08
150  ST+ 10
151  ST+ 11
152  XEQ 00
153  LBL 06
154  RCL IND 07
155  RCL 02
156  *		 157  STO IND 12
158  RCL 20
159  *
160  RCL IND 11
161  RCL 21
162  *
163  +
164  RCL IND 10
165  RCL 22
166  *
167  +
168  RCL IND 08
169  RCL 23
170  *
171  +
172  ST+ IND 06
173  DSE 12
174  DSE 11
175  DSE 10
176  DSE 08
177  DSE 07
178  DSE 06
179  GTO 06
180  RCL 01
181  ST+ 06
182  ST+ 07
183  ST+ 08
184  ST+ 10
185  ST+ 11
186  ST+ 12
187  RCL 02
188  2
189  /
190  RCL 04
191  +
192  STO 50
193  XEQ 00
194  LBL 07
195  RCL IND 07
196  RCL 02
197  *
198  STO IND 13
199  RCL 24
200  *
201  RCL IND 12
202  RCL 25
203  *
204  +
205  RCL IND 11
206  RCL 26
207  *
208  +
209  RCL IND 10
210  RCL 27
211  *
212  +
213  RCL IND 08
214  RCL 28
215  *
216  +
217  ST+ IND 06
218  DSE 13
219  DSE 12
220  DSE 11
221  DSE 10
222  DSE 08
223  DSE 07
224  DSE 06
225  GTO 07
226  RCL 01
227  ST+ 06
228  ST+ 07
229  ST+ 08
230  ST+ 10
231  ST+ 11
232  ST+ 12
233  ST+ 13
234  RCL 02		 235  RCL 29
236  *
237  RCL 04
238  +
239  STO 50
240  XEQ 00
241  LBL 08
242  RCL IND 07
243  RCL 02
244  *
245  STO IND 09
246  9
247  /
248  RCL IND 13
249  RCL 30
250  *
251  +
252  RCL IND 12
253  RCL 31
254  *
255  +
256  RCL IND 08
257  14
258  /
259  +
260  ST+ IND 06
261  DSE 13
262  DSE 12
263  DSE 09
264  DSE 08
265  DSE 07
266  DSE 06
267  GTO 08
268  RCL 01
269  ST+ 06
270  ST+ 07
271  ST+ 08
272  ST+ 09
273  ST+ 12
274  ST+ 13
275  XEQ 00
276  LBL 09
277  RCL IND 07
278  RCL 02
279  *
280  STO IND 10
281  RCL 32
282  *
283  RCL IND 09
284  RCL 33
285  *
286  +
287  RCL IND 13
288  RCL 34
289  *
290  +
291  RCL IND 12
292  RCL 35
293  *
294  +
295  RCL IND 08
296  32
297  /
298  +
299  ST+ IND 06
300  DSE 13
301  DSE 12
302  DSE 10
303  DSE 09
304  DSE 08
305  DSE 07
306  DSE 06
307  GTO 09
308  RCL 01
309  ST+ 06
310  ST+ 07
311  ST+ 08
312  ST+ 09		 313  ST+ 10
314  ST+ 12
315  ST+ 13
316  RCL 02
317  2
318  /
319  RCL 04
320  +
321  STO 50
322  XEQ 00
323  LBL 10
324  RCL IND 07
325  RCL 02
326  *
327  STO IND 11
328  RCL 36
329  *
330  RCL IND 10
331  RCL 37
332  *
333  +
334  RCL IND 09
335  RCL 38
336  *
337  +
338  RCL IND 13
339  RCL 39
340  *
341  +
342  RCL IND 12
343  9
344  /
345  +
346  RCL IND 08
347  14
348  /
349  +
350  ST+ IND 06
351  DSE 13
352  DSE 12
353  DSE 11
354  DSE 10
355  DSE 09
356  DSE 08
357  DSE 07
358  DSE 06
359  GTO 10
360  RCL 01
361  ST+ 06
362  ST+ 07
363  ST+ 08
364  ST+ 09
365  ST+ 10
366  ST+ 11
367  ST+ 12
368  ST+ 13
369  RCL 02
370  RCL 16
371  *
372  RCL 04
373  +
374  STO 50
375  XEQ 00
376  LBL 11
377  RCL IND 07
378  RCL 02
379  *
380  X<> IND 12
381  RCL 40
382  *
383  RCL IND 12
384  RCL 41
385  *
386  +
387  RCL IND 11
388  RCL 42
389  *
390  +		 391  RCL IND 10
392  RCL 43
393  *
394  +
395  RCL IND 09
396  RCL 44
397  *
398  +
399  RCL IND 13
400  RCL 45
401  *
402  +
403  ST+ IND 06
404  DSE 13
405  DSE 12
406  DSE 11
407  DSE 10
408  DSE 09
409  DSE 07
410  DSE 06
411  GTO 11
412  RCL 01
413  ST+ 06
414  ST+ 07
415  ST+ 09
416  ST+ 10
417  ST+ 11
418  ST+ 12
419  ST+ 13
420  RCL 02
421  RCL 04
422  +
423  STO 50
424  XEQ 00
425  LBL 12
426  RCL IND 07
427  RCL 02
428  *
429  RCL IND 08
430  +
431  RCL 46
432  *
433  RCL IND 10
434  RCL IND 12
435  +
436  RCL 47
437  *
438  +
439  RCL IND 11
440  RCL 48
441  *
442  +
443  ST+ IND 06
444  DSE 12
445  DSE 11
446  DSE 10
447  DSE 08
448  DSE 07
449  DSE 06
450  GTO 12
451  RCL 01
452  ST+ 06
453  ST+ 07
454  ST+ 08
455  ST+ 10
456  ST+ 11
457  ST+ 12
458  DSE 49
459  GTO 01 
460  RCL 53
461  RCL 52
462  RCL 51
463  RCL 50
464  END
 

 
   ( 765 bytes / SIZE 51+9n )
 
 

      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
 

  LBL "T"      *                 RCL 51      -                   LASTX      RCL 53       RTN            "T"  ASTO 00
  RCL 51      *                 RCL 52      RCL 50        RCL 51      *
  RCL 52      CHS           +                 *                  *                 -
  RCL 53      STO 54      RCL 53       STO 55        RCL 52      STO 56

 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 "RK8"  >>>>      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 )
-If you replace lines 431 thru 442 by  9  *  RCL IND 10  RCL IND 12  +  49  *  +  RCL IND 11  64  *  +  PI  R-D  /
  registers  R46-R47-R48  will be unused.
 
 

     d) Runge-Kutta-Fehlberg 4-5 Method
 

-The following program uses a 4th-order formula to compute the solution,
 and the difference between this 4th-order formula and a 5th-order formula provides an error estimate.
 

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

                                      •  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

-If you don't have an HP-41 CX, replace line 31 by   0   LBL 00   STO IND Y   ISG Y   GTO 00
-Lines 290 & 298 are three-byte GTOs
 
 

  01  LBL "RKF" 
  02  21.021
  03  RCL 01
  04   E6
  05  /
  06  +
  07  RCL 01
  08  7
  09  *
  10  +
  11  STO 05
  12  FRC
  13  RCL 01
  14  .1
  15  %
  16  +
  17  +
  18  RCL 01
  19  +
  20  21
  21  +
  22  STO 06
  23  RCL 01
  24  ST+ X
  25   E3
  26  /
  27  RCL 01
  28  +
  29  21.02
  30  +
  31  CLRGX
  32  RCL 01
  33  RCL 01
  34  LASTX
  35  1
  36  -
  37  +
  38  STO 07
  39  +
  40  STO 08
  41  +
  42  STO 13
  43  +
  44  STO 09
  45  +
  46  STO 10
  47  +
  48  STO 11
  49  +
  50  STO 12		   51  LBL 10
  52  RCL 03
  53  STO 04
  54  RCL 06
  55  REGSWAP
  56  LBL 01
  57  RCL 20
  58  STO 14
  59  RCL 05
  60  REGSWAP
  61  REGMOVE
  62  XEQ IND 00
  63  LBL 02
  64  RCL IND 08
  65  RCL 02
  66  *
  67  STO IND 09
  68  ST+ X
  69  9
  70  /
  71  ST+ IND 07
  72  DSE 09
  73  DSE 08
  74  DSE 07
  75  GTO 02
  76  RCL 01
  77  ST+ 07
  78  ST+ 08
  79  ST+ 09
  80  RCL 02
  81  ST+ X
  82  LASTX
  83  /
  84  ST+ 20
  85  XEQ IND 00
  86  RCL 05
  87  REGMOVE
  88  LBL 03
  89  RCL IND 08
  90  RCL 02
  91  *
  92  4
  93  /
  94  STO IND 10
  95  RCL IND 09
  96  12
  97  /
  98  +
  99  ST+ IND 07
100  DSE 10		 101  DSE 09
102  DSE 08
103  DSE 07
104  GTO 03
105  RCL 01
106  ST+ 07
107  ST+ 08
108  ST+ 09
109  ST+ 10
110  RCL 02
111  9
112  /
113  ST+ 20
114  XEQ IND 00
115  RCL 05
116  REGMOVE
117  LBL 04
118  RCL IND 08
119  RCL 02
120  *
121  STO IND 11
122  270
123  *
124  RCL IND 10
125  972
126  *
127  -
128  RCL IND 09
129  69
130  *
131  +
132  128
133  /
134  ST+ IND 07
135  DSE 11
136  DSE 10
137  DSE 09
138  DSE 08
139  DSE 07
140  GTO 04
141  RCL 01
142  ST+ 07
143  ST+ 08
144  ST+ 09
145  ST+ 10
146  ST+ 11
147  RCL 02
148  2.4
149  /
150  ST+ 20		 151  XEQ IND 00
152  RCL 05
153  REGMOVE
154  LBL 05
155  RCL IND 08
156  RCL 02
157  *
158  64
159  *
160  STO IND 12
161  RCL IND 09
162  85
163  *
164  -
165  60
166  /
167  RCL IND 10
168  RCL IND 11
169  5
170  /
171  -
172  27
173  *
174  +
175  ST+ IND 07
176  DSE 12
177  DSE 11
178  DSE 10
179  DSE 09
180  DSE 08
181  DSE 07
182  GTO 05
183  RCL 01
184  ST+ 07
185  ST+ 08
186  ST+ 09
187  ST+ 10
188  ST+ 11
189  ST+ 12
190  RCL 02
191  ST+ 14
192  RCL 14
193  STO 20
194  XEQ IND 00
195  RCL 05
196  REGMOVE
197  LBL 06
198  RCL IND 08
199  RCL 02
200  *		 201  15
202  *
203  ST+ IND 12
204  RCL IND 12
205  RCL IND 11
206  351
207  *
208  +
209  RCL IND 10
210  540
211  *
212  -
213  RCL IND 09
214  65
215  *
216  +
217  432
218  /
219  ST+ IND 07
220  DSE 12
221  DSE 11
222  DSE 10
223  DSE 09
224  DSE 08
225  DSE 07
226  GTO 06
227  RCL 01
228  ST+ 07
229  ST+ 08
230  ST+ 09
231  ST+ 10
232  ST+ 11
233  ST+ 12
234  RCL 02
235  6
236  /
237  ST- 20
238  XEQ IND 00
239  RCL 05
240  REGMOVE
241  LBL 07
242  RCL IND 08
243  RCL 02
244  *
245  72
246  *
247  RCL IND 12
248  -
249  RCL IND 11
250  9		 251  *
252  +
253  RCL IND 09
254  ST+ X
255  -
256  3
257  1/X
258  %
259  ST- IND 13
260  RCL IND 12
261  RCL IND 11
262  81
263  *
264  +
265  PI
266  R-D
267  /
268  RCL IND 09
269  9
270  /
271  +
272  ST+ IND 07
273  DSE 13
274  DSE 12
275  DSE 11
276  DSE 09
277  DSE 08
278  DSE 07
279  GTO 07
280  RCL 01
281  ST+ 07
282  ST+ 08
283  ST+ 09
284  ST+ 11
285  ST+ 12
286  ST+ 13
287  RCL 14
288  STO 20
289  DSE 04
290  GTO 01
291  RCL 06
292  REGSWAP
293  RCL 23
294  RCL 22
295  RCL 21
296  RCL 20
297  RTN
298  GTO 10
299  END

 
   ( 481 bytes / SIZE 21+8n )
 
 

      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
 

  LBL "T"      *                 RCL 21      -                   RCL 20      RCL 23       RTN            "T"  ASTO 00
  RCL 21      *                 RCL 22      RCL 20        RCL 21      *
  RCL 22      CHS           +                 *                  *                 -
  RCL 23      STO 24      RCL 23       STO 25        RCL 22      STO 26

 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 )

-If you want to use the 5th-order formula to compute yi(x) , replace line 259 with ST+ IND 07   ST+ IND 13  but the errors will be probably overestimated.
-In the above example, it gives:

     y(1) = 0.258207897
     z(1) = 1.157624056
     u(1) = 0.842178340
 

     e)  Bulirsch-Stoer Method
 

-This program uses the formulae given in "The Bulirsch-Stoer Method for the HP-41" and - first of all - in "Numerical Recipes"
 

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

                                      •  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

-Lines 194-245-252-262-266  are three-byte  GTOs
 
 

  01  LBL "BST"
  02  STO 03
  03  1
  04  STO 04
  05  9
  06  STO 05
  07  21.021
  08  RCL 01
  09   E6
  10  /
  11  +
  12  RCL 01
  13  5
  14  *
  15  +
  16  STO 08
  17  RCL 01
  18  ST+ X
  19  ST- Y
  20  .1
  21  %
  22  -
  23  -
  24  STO 16       
  25  20.02
  26  RCL 01
  27  +
  28  STO 10
  29  INT
  30  .1
  31  %
  32  RCL 01
  33  STO T
  34  +
  35  +
  36  STO 11
  37  +
  38  STO 12
  39  +
  40  STO 13
  41  +
  42  STO 14
  43  GTO 01
  44  LBL 00
  45  2		   46  ST/ 04
  47  GTO 02
  48  LBL 01
  49  RCL 20
  50  STO 15
  51  RCL 08
  52  REGSWAP
  53  REGMOVE
  54  XEQ IND 00
  55  RCL 16
  56  REGMOVE
  57  6
  58  RCL 05
  59  X>Y?
  60  GTO 02
  61  RCL 04
  62  ST+ 04
  63  LBL 02
  64  SF 09
  65  RCL 04
  66  STO 17
  67  ABS
  68  RCL 03
  69  RCL 15
  70  -
  71  ABS
  72  X<=Y?
  73  CF 09
  74  X>Y?
  75  X<>Y
  76  LASTX
  77  SIGN
  78  *
  79  STO 04
  80  SF 10
  81  CLX
  82  STO 05
  83  RCL 08
  84  INT
  85  STO 09
  86  LBL 03
  87  RCL 08
  88  REGMOVE 
  89  2
  90  ST+ 05		   91  18
  92  RCL 05
  93  X=Y?
  94  GTO 00
  95  RCL 01
  96  ST+ 09
  97  SIGN
  98  -
  99   E3
100  /
101  ISG X
102  STO 07
103  RCL 04
104  RCL 05
105  /
106  STO 06
107  LBL 04
108  RCL IND 12
109  RCL 06
110  *
111  RCL IND 10
112  STO IND 13
113  +
114  STO IND 14
115  DSE 14
116  DSE 13
117  DSE 12
118  DSE 10
119  GTO 04
120  RCL 01
121  ST+ 10
122  ST+ 12
123  ST+ 13
124  ST+ 14
125  LBL 05
126  RCL 08
127  RCL 01
128  -
129  REGMOVE
130  RCL 06
131  RCL07
132  INT
133  *
134  RCL 15
135  +		 136  STO 20
137  XEQ IND 00
138  LBL 06
139  RCL IND 11
140  RCL 06
141  ST+ X
142  *
143  RCL IND 14
144  X<> IND 13
145  +
146  STO IND 14
147  DSE 14
148  DSE 13
149  DSE 11
150  GTO 06
151  RCL 01
152  ST+ 11
153  ST+ 13
154  ST+ 14
155  ISG 07
156  GTO 05
157  RCL 08
158  RCL 01
159  -
160  REGMOVE
161  RCL 04
162  RCL 15
163  +
164  STO 20
165  XEQ IND 00
166  LBL 07
167  RCL IND 11
168  RCL 06
169  *
170  RCL IND 13
171  +
172  ST+ IND 10
173  2
174  ST/ IND 10
175  SIGN
176  ST- 11
177  ST- 13
178  DSE 10
179  GTO 07
180  21		 181  RCL 09
182   E3
183  /
184  +
185  RCL 01
186  ST+ 10
187  ST+ 11
188  ST+ 13
189   E6
190  /
191  +
192  REGMOVE
193  FS?C 10
194  GTO 03  
195  RCL 09
196  RCL 01
197  +
198  DSE X
199   E3
200  /
201  RCL 09
202  +
203  STO 06
204  CLX
205  STO 18
206  LBL 08
207  RCL 05
208  2
209  ST- Y
210  /
211  STO 07
212  LBL 09
213  RCL 06
214  RCL 07
215  RCL 01
216  *
217  -
218  RCL IND 06
219  ENTER^
220  X<> IND Z
221  STO Z
222  -
223  RCL 05
224  X^2
225  ST* Y		 226  RCL 07
227  ST+ X
228  X^2
229  -
230  /
231  +
232  STO IND 06
233  DSE 07
234  GTO 09
235  LASTX
236  ABS
237  RCL 18
238  X<Y?
239  X<>Y
240  STO 18
241  ISG 06
242  GTO 08
243  RCL 02
244  X<Y?
245  GTO 03  
246  RCL 09
247  RCL 08
248  FRC
249  +
250  REGMOVE 
251  FS? 09
252  GTO 01  
253  RCL 23
254  RCL 22
255  RCL 21
256  RCL 20
257  RTN
258  STO 03
259  RCL 04
260  RCL 17
261  X=Y?
262  GTO 01 
263  STO 04
264  9
265  STO 05
266  GTO 01
267  END

 
   ( 395 bytes / SIZE 21+14n )
 
 

      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
 

  LBL "T"      *                 RCL 21      -                   LASTX      RCL 23       RTN            "T"  ASTO 00
  RCL 21      *                 RCL 22      RCL 20        RCL 21      *
  RCL 22      CHS           +                 *                  *                 -
  RCL 23      STO 24      RCL 23       STO 25        RCL 22      STO 26

 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 "BST"  >>>>    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 )
-If you want to choose your own initial stepsize, place it in Y-register after replacing line 03 by X<>Y
-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.
-If this is not satisfied, even with h/16, the stepsize is divided by 2 ( LBL 00 )
-If - for instance - you want to stop the extrapolation with  h/14, replace line 91 by 16 ( instead of 18 )
-Thus, the SIZE may be reduced to 21+13n
-On the other hand, if you want to continue the extrapolation with  h/18, replace line 91 by 20
-In this case, the SIZE must be at least 21+15n
 

2°) Systems of second-order Conservative Equations
 

-The 3 following programs use the 4th-order Runge-Kutta formula:

   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)
 

     a)  Systems of 2 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
 

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

                               •   R01 = x0      •   R04 = h/2 = half of the stepsize       •   R06 = y'0          R08 to R12: 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.

-Line 96 is a three-byte  GTO 01
 
 

  01  LBL "RKS2"
  02  RCL 05
  03  STO 12
  04  LBL 01
  05  RCL 03
  06  RCL 02
  07  RCL 01
  08  XEQ IND 00
  09  RCL 04
  10  ST+ 01
  11  ST* Z
  12  *
  13  STO 09
  14  2
  15  /
  16  RCL 07
  17  +
  18  RCL 04
  19  *
  20  RCL 03		   21  +
  22  X<>Y
  23  STO 08
  24  2
  25  /
  26  RCL 06
  27  +
  28  RCL 04
  29  *
  30  RCL 02
  31  +
  32  RCL 01
  33  XEQ IND 00
  34  RCL 04
  35  ST+ 01
  36  ST* Z
  37  *
  38  STO 11
  39  RCL 07
  40  +		   41  X<>Y
  42  STO 10
  43  RCL 06
  44  +
  45  RCL 04
  46  ST+ X
  47  ST* Z
  48  *
  49  RCL 03
  50  ST+ Z
  51  CLX
  52  RCL 02
  53  +
  54  RCL 01
  55  XEQ IND 00
  56  RCL 04
  57  ST* Z
  58  *
  59  RCL 09
  60  +		   61  RCL 11
  62  ST+ X
  63  ST+ 09
  64  ST+ X
  65  +
  66  3
  67  ST/ 09
  68  /
  69  X<> 07
  70  ST+ 07
  71  RCL 09       
  72  +
  73  X<>Y
  74  RCL 08
  75  +
  76  RCL 10
  77  ST+ X
  78  ST+ 08
  79  ST+ X
  80  +		   81  3
  82  ST/ 08
  83  /
  84  X<> 06
  85  ST+ 06
  86  RCL 08        
  87  +
  88  RCL 04
  89  ST+ X
  90  ST* Z
  91  *
  92  ST+ 02
  93  X<>Y
  94  ST+ 03
  95  DSE 12
  96  GTO 01
  97  RCL 03
  98  RCL 02
  99  RCL 01
100  END

 
   ( 139 bytes / SIZE 013 )
 
 

      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"      CHS        *                     "T" ASTO 00
   RDN          ST* Z       RTN
   STO Z        -
   X<>Y        R^

    0     STO 01              1  STO 06
    2     STO 02              1  STO 07
    1     STO 03
  0.05  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 )
 

     b)  Systems of 3 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/2 = half of the stepsize       •   R07 = y'0
                               •   R02 = y0      •   R06 = N = number of steps            •   R08 = z'0
                               •   R03 = z0                   R10 to R16: 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.

-Line 136 is a three-byte  GTO 01
 
 

  01  LBL "RKS3"
  02  RCL 06
  03  STO 16
  04  LBL 01
  05  RCL 04
  06  RCL 03
  07  RCL 02
  08  RCL 01
  09  XEQ IND 00
  10  RCL 05
  11  ST+ 01
  12  ST* Z
  13  ST* T
  14  *
  15  STO 12
  16  2
  17  /
  18  RCL 09
  19  +
  20  RCL 05
  21  *
  22  RCL 04
  23  +
  24  X<>Y
  25  STO 11
  26  2
  27  /
  28  RCL 08
  29  +		   30  RCL 05
  31  *
  32  RCL 03
  33  +
  34  R^
  35  STO 10
  36  2
  37  /
  38  RCL 07
  39  +
  40  RCL 05
  41  *
  42  RCL 02
  43  +
  44  RCL 01
  45  XEQ IND 00
  46  RCL 05
  47  ST+ 01
  48  ST* Z
  49  ST* T
  50  *
  51  STO 15
  52  RCL 09
  53  +
  54  X<>Y
  55  STO 14
  56  RCL 08
  57  +
  58  R^		   59  STO 13
  60  RCL 07
  61  +
  62  RCL 05
  63  ST+ X
  64  ST* Z
  65  ST* T
  66  *
  67  RCL 04
  68  ST+ T
  69  CLX
  70  RCL 03
  71  ST+ Z
  72  CLX
  73  RCL 02
  74  +
  75  RCL 01
  76  XEQ IND 00
  77  RCL 05
  78  ST* Z
  79  ST* T
  80  *
  81  RCL 12
  82  +
  83  RCL 15
  84  ST+ X
  85  ST+ 12
  86  ST+ X
  87  +		   88  3
  89  ST/ 12
  90  /
  91  X<> 09
  92  ST+ 09
  93  RCL 12
  94  +
  95  X<>Y
  96  RCL 11
  97  +
  98  RCL 14
  99  ST+ X
100  ST+ 11
101  ST+ X
102  +
103  3
104  ST/ 11
105  /
106  X<> 08        
107  ST+ 08
108  RCL 11
109  +
110  R^
111  RCL 10
112  +
113  RCL 13
114  ST+ X
115  ST+ 10
116  ST+ X		 117  +
118  3
119  ST/ 10
120  /
121  X<> 07
122  ST+ 07
123  RCL 10
124  +
125  RCL 05
126  ST+ X
127  ST* Z
128  ST* T
129  *
130  ST+ 02
131  X<>Y
132  ST+ 03
133  R^
134  ST+ 04
135  DSE 16        
136  GTO 01 
137  RCL 04
138  RCL 03
139  RCL 02
140  RCL 01
141  END
 

 
   ( 194 bytes / SIZE 017 )
 
 

      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^               RCL Z         X<> Z         LASTX                 "T"  ASTO 00
  R^               ST* Y         R^                ST+ Y         *
  CHS            ST+ L         ST+ L          ST* Z          X<>Y
  STO L         CLX           ST* Y          X<> T         RTN

  0     STO 01
  1     STO 02  STO 03  STO 07  STO 08  STO 09
  2     STO 04
0.05  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 "RKS" below, for 3 equations.
 

     c)  Systems of n Equations
 
 

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

                                      •  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

-Line 140 is a three-byte  GTO 01
 
 

  01  LBL "RKS"
  02  RCL 03
  03  STO 04
  04  11.011
  05  RCL 01
  06   E6
  07  /
  08  +
  09  RCL 01
  10  5
  11  *
  12  +
  13  STO 05
  14  RCL 01
  15  RCL 01
  16  RCL 01
  17  10.01
  18  +
  19  STO 06
  20  +
  21  STO 07
  22  +
  23  STO 08
  24  +
  25  STO 09
  26  LBL 01
  27  RCL 05
  28  REGSWAP
  29  REGMOVE 		   30  FRC
  31  11
  32  +
  33  RCL 01
  34  .1
  35  %
  36  ST+ X
  37  +
  38  +
  39  REGMOVE
  40  XEQ IND 00
  41  LBL 02
  42  RCL IND 07
  43  RCL 02
  44  *
  45  STO IND 09
  46  4
  47  /
  48  RCL IND 08
  49  +
  50  RCL 02
  51  *
  52  2
  53  /
  54  ST+ IND 06
  55  DSE 09
  56  DSE 08
  57  DSE 07
  58  DSE 06		   59  GTO 02
  60  RCL 01
  61  ST+ 06
  62  ST+ 07
  63  ST+ 08
  64  ST+ X
  65  ST+ 09
  66  RCL 02
  67  LASTX
  68  /
  69  ST+ 10
  70  XEQ IND 00
  71  RCL 05
  72  REGMOVE
  73  LBL 03
  74  RCL IND 07
  75  RCL 02
  76  *
  77  STO IND 09
  78  2
  79  /
  80  RCL IND 08
  81  +
  82  RCL 02
  83  *
  84  ST+ IND 06
  85  DSE 09
  86  DSE 08
  87  DSE 07		   88  DSE 06
  89  GTO 03
  90  RCL 01
  91  ST+ 06
  92  ST+ 08
  93  ST+ 09
  94  3
  95  *
  96  ST+ 07
  97  RCL 02
  98  2
  99  /
100  ST+ 10
101  XEQ IND 00
102  RCL 05
103  REGMOVE
104  LBL 04
105  RCL IND 09
106  ST+ X
107  ST+ IND 07
108  X<> IND 07
109  ST+ IND 07
110  6
111  /
112  RCL IND 08
113  +
114  RCL 02
115  *
116  ST+ IND 06		 117  RCL 08
118  RCL 01
119  -
120  RCL IND 07
121  RCL IND Y
122  RCL 02
123  *
124  +
125  6
126  /
127  RCL IND 08
128  +
129  STO IND Y
130  DSE 09
131  DSE 08
132  DSE 07
133  DSE 06
134  GTO 04
135  RCL 01
136  ST+ 06
137  ST- 07
138  ST+ 08
139  DSE 04
140  GTO 01 
141  RCL 13
142  RCL 12
143  RCL 11
144  RCL 10
145  END

 
   ( 225 bytes / SIZE 6n+11 )
 
 

      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
 

  LBL "T"      *                 RCL 11      -                   RCL 10      RCL 13       RTN            "T"  ASTO 00
  RCL 11      *                 RCL 12      RCL 10        RCL 11      *
  RCL 12      CHS           +                 *                  *                 -
  RCL 13      STO 14      RCL 13       STO 15        RCL 12      STO 16

 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 "RKS"  >>>>    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 needed )
 

References:

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