The Museum of HP Calculators

# Systems of Differential Equations for the HP-41

This program is Copyright © 2007 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

## Overview

1°) Systems of first-order Differential Equations

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

2°) Systems of second-order Conservative Equations

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

-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  = 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

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 gives:            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

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                  a three-byte GTO  ( or replace this line by  GTO 15  and line 32 by  LBL 15 )
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

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                  a three-byte GTO
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

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                  If you don't have an HP-41 CX, replace this line by   0   LBL 00   STO IND Y   ISG Y   GTO 00
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                  a three-byte GTO
291  RCL 06
292  REGSWAP
293  RCL 23
294  RCL 22
295  RCL 21
296  RCL 20
297  RTN
298  GTO 10                  a three-byte GTO
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

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                  a three-byte GTO
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                  a three-byte GTO
246  RCL 09
247  RCL 08
248  FRC
249  +
250  REGMOVE
251  FS? 09
252  GTO 01                  a three-byte GTO
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                  a three-byte GTO
263  STO 04
264  9
265  STO 05
266  GTO 01                  a three-byte GTO
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, rounded to 9D are          z(2) = 3.886706159
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.

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                  a three-byte GTO
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          R10 to R16: temp
•   R02 = y0      •   R06 = N = number of steps            •   R08 = z'0
•   R03 = z0                                                                 •   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.

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                  a three-byte GTO
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

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                  a three-byte GTO
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