The Museum of HP Calculators

# More Implicit-Runge-Kutta methods

This program is 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

-Two kinds of programs are presented below:

1°) IRK is a general Runge-Kutta program suitable to all implicit formulae
2°) I2RK is slightly less general.

-These methods are applied to solve:

a) A first order differential equation:       dy/dx = f(x,y)                                          with the initial value:   y(x0) = y0
b) A system of 2 differential equations:   dy/dx = f(x,y,z)  ;   dz/dx = g(x,y,z)          -------------------    y(x0) = y0  ;   z(x0) = z0

1°)  General implicit Runge-Kutta methods

A) First order differential equations:  dy/dx = f (x;y)

-A n-stage implicit Runge-Kutta method is defined by   n(n+2)  coefficients ai,j ; bi ; ci

k1 = h.f ( x + b1h ; y + a1,1 k1 +  a1,2 k2 +  ..... + a1,n kn )
.....................................................................................            (  actually,     ai,1 + ai,2 + .......... + ai,n = bi  )
kn = h.f ( x + bnh ; y + an,1 k1 +  an,2 k2 +  ..... + an,n kn )

then:     y(x+h) = y(x) + c1.k1+ ................ + cn.kn

-The following program will work for any implicit method.

Data Registers:     Registers R00 to R05 and    Rn+11 to Rn2+3n+10     are to be initialized before executing "IRK"

•   R00:  f name   •   R01 = x0                                      R06 to R10:  temp             •  Rn+11 to Rn2+3n+10 =  the ( n2 + 2n ) coefficients ai,j ; bi ; ci
•   R02 = y0                                      R11 = k1                                                                     which are to be stored as shown below:
•   R03 = h = stepsize                        R12 = k2
•   R04 = N = number of steps          ............
•   R05 = n  = number of stages        Rn+10 = kn

•  Rn+11   = a1,1        •  Rn+12   = a1,2   ...........................    •  R2n+10 = a1,n     •  R2n+11 = b1
•  R2n+12 = a2,1        •  R2n+13 = a2,2   ...........................    •  R3n+11 = a2,n     •  R3n+12 = b2

..................................................................................................................................................

•  Rn2+n+10 = an,1     •  Rn2+n+11 = an,2   .......................   •  Rn2+2n+9 = an,n   •  Rn2+2n+10 = bn

•  Rn2+2n+11 = c1      •  Rn2+2n+12 = c2   .......................   •  Rn2+3n+10 = cn

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

>>>>>       In other words, this subroutine has to change the stack from   Y = y
X = x        into       X' = f(x;y)

01  LBL "IRK"
02  RCL 04
03  STO 06
04  RCL 05
05  10
06  +
07  .1
08  %
09  11
10  +
11  STO 07
12  CLRGX               if you don't have an HP-41CX , replace line 12 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00
13  LBL 01
14  RCL 05
15  STO 08
16  RCL 07
17  FRC
18  11.9
19  ST+ 08
20  INT
21  +
22  STO 07
23  CLX
24  STO 10
25  LBL 02
26  RCL 07
27  FRC
28  11
29  +
30  STO 09
31  CLX
32  LBL 03
33  RCL IND 08
34  RCL IND 09
35  *
36  +
37  ISG 08
38  ISG 09
39  GTO 03
40  RCL 02
41  +
42  RCL 03
43  RCL IND 08
44  *
45  RCL 01
46  +
47  XEQ IND 00
48  RCL 03
49  *
50  ENTER^
51  X<> IND 07
52  -
53  ABS
54  ST+ 10
55  ISG 08
56  ISG 07
57  GTO 02
58  RCL 10
59  VIEW X
60   E-9                     or another "small" number
61  X<Y?
62  GTO 01
63  RCL 05
64  ST- 07
65  CLX
66  LBL 04
67  RCL IND 07
68  RCL IND 08
69  *
70  +
71  ISG 08
72  ISG 07
73  GTO 04
74  ST+ 02
75  RCL 03
76  ST+ 01
77  DSE 06
78  GTO 01
79  RCL 02
80  RCL 01
81  CLD
82  END

( 126 bytes / SIZE n2+3n+11 )

 STACK INPUTS OUTPUTS Y / y(x0+N.h) X / x0+N.h

-Implicit methods are very accurate:  n-stage 2n-order methods do exist which would be impossible with explicit methods,
but, on the other hand, we have to solve a nxn non-linear system within each step!
-For instance, here are the 48 coefficients of a 6-stage 12-order method, based on Gauss-Legendre quadrature.
-These coefficients are given with 20 decimals in case you want to use them with a calculator working with more than 10 digits:

R17 = a1,1 =  0.04283 11230 94792 58626          R24 = a2,1 =  0.09267 34914 30378 86319          R31 = a3,1 =  0.08224 79226 12843 87381
R18 = a1,2 = -0.01476 37259 97197 41248          R25 = a2,2 =  0.09019 03932 62034 65189          R32 = a3,2 =  0.19603 21623 33245 00606
R19 = a1,3 =  0.00932 50507 06477 75119          R26 = a2,3 = -0.02030 01022 93239 58595          R33 = a3,3 =  0.11697 84836 43172 76185
R20 = a1,4 = -0.00566 88580 49483 51191          R27 = a2,4 =  0.01036 31562 40246 42374          R34 = a3,4 = -0.02048 25277 45656 09763
R21 = a1,5 =  0.00285 44333 15099 33514          R28 = a2,5 = -0.00488 71929 28037 67146          R35 = a3,5 =  0.00798 99918 99662 33580
R22 = a1,6 = -0.00081 27801 71264 76211          R29 = a2,6 =  0.00135 55610 55485 06178          R36 = a3,6 = -0.00207 56257 84866 33419
R23 =  b1  =  0.03376 52428 98423 98609          R30 =  b2  =  0.16939 53067 66867 74317          R37 =   b3 =   0.38069 04069 58401 54568

R38 = a4,1 =  0.08773 78719 74451 50671          R45 = a5,1 =  0.08430 66851 34100 11074          R52 = a6,1 =  0.08647 50263 60849 93463
R39 = a4,2 =  0.17239 07946 24406 96799          R46 = a5,2 =  0.18526 79794 52106 97525          R53 = a6,2 =  0.17752 63532 08969 96865
R40 = a4,3 =  0.25443 94950 32001 62132          R47 = a5,3 =  0.22359 38110 46099 09996          R54 = a6,3 =  0.23962 58253 35829 03560
R41 = a4,4 =  0.11697 84836 43172 76185          R48 = a5,4 =  0.25425 70695 79585 10965          R55 = a6,4 =  0.22463 19165 79867 77250
R42 = a4,5 = -0.01565 13758 09175 70227         R49 = a5,5 =  0.09019 03932 62034 65189           R56 = a6,5 =  0.19514 45125 21266 71626
R43 = a4,6 =  0.00341 43235 76741 29871          R50 = a5,6 = -0.00701 12452 40793 69067          R57 = a6,6 =  0.04283 11230 94792 58626
R44 =  b4  =  0.61930 95930 41598 45432          R51 =  b5  =  0.83060 46932 33132 25683          R58 =  b6  =  0.96623 47571 01576 01391

R59 =  c1  =  0.08566 22461 89585 17252          R61 =  c3  =  0.23395 69672 86345 52369           R63 =  c5  =  0.18038 07865 24069 30378
R60 =  c2  =  0.18038 07865 24069 30378          R62 =  c4  =  0.23395 69672 86345 52369           R64 =  c6  =  0.08566 22461 89585 17252

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

-Initialize registers R17 thru R64

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.5
0.5  STO 03
2    STO 04    ( the number of steps )
6    STO 05    ( we are using a 6-stage method )                        XEQ "IRK"   yields          x = 1                         ( in 5mn36s )
X<>Y                  y(1) = 0.3678794411     ( error = -10-10 )
( With  h = N = 1 the error is only  -8.4 10-9 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | kn - k'n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

B) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

-The same formulae may be generalized to a system of 2 differential equations, and we have to solve a 2nx2n non-linear system within each step.

Data Registers:     Registers R00 to R06 and    R2n+14 to Rn2+4n+13     are to be initialized before executing "IRKB"

•   R00:  f name   •   R01 = x0                                      R07 to R13:  temp             •  R2n+14 to Rn2+4n+13 =  the ( n2 + 2n ) coefficients ai,j ; bi ; ci
•   R02 = y0                                      R14 to R2n+13 = ki                                                      which are to be stored as shown below:
•   R03 = z0
•   R04 = h   = stepsize
•   R05 = N  = number of steps
•   R06 = n   = number of stages

•  R2n+14 = a1,1        •  R2n+15 = a1,2   ...........................    •  R3n+13 = a1,n         •  R3n+14 = b1
•  R3n+15 = a2,1        •  R3n+16 = a2,2   ...........................    •  R4n+14 = a2,n         •  R4n+15 = b2

.....................................................................................................................................................

•  Rn2+2n+13 = an,1   •  Rn2+2n+14 = an,2   ......................   •  Rn2+3n+12 = an,n   •  Rn2+3n+13 = bn

•  Rn2+3n+14 = c1      •  Rn2+3n+15 = c2   ........................   •  Rn2+4n+13 = cn

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

Z = z
>>>>>    In other words, this subroutine has to change the stack from   Y = y       into     Y' = f(x;y;z)
X = x                  X' = g(x;y;z)

001  LBL "IRKB"
002  RCL 05
003  STO 07
004  RCL 06
005  ST+ X
006  13
007  +
008  .1
009  %
010  14
011  +
012  STO 09
013  CLRGX               if you don't have an HP-41CX , replace line 13 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00
014  LBL 01
015  RCL 09
016  FRC
017  14.9
018  STO 10
019  INT
020  +
021  STO 08
022  RCL 06
023  ST+ 10
024  ST+ 10
025  +
026  STO 09
027  CLX
028  STO 13
029  LBL 02
030  RCL 09
031  FRC
032  14
033  +
034  STO 11
035  RCL 06
036  +
037  STO 12
038  CLST
039  LBL 03
040  X<>Y
041  RCL IND 10
042  RCL IND 11
043  *
044  +
045  X<>Y
046  RCL IND 10
047  RCL IND 12
048  *
049  +
050  ISG 10
051  ISG 11
052  ISG 12
053  GTO 03
054  RCL 03
055  +
056  X<>Y
057  RCL 02
058  +
059  RCL 04
060  RCL IND 10
061  *
062  RCL 01
063  +
064  XEQ IND 00
065  RCL 04
066  ST* Z
067  *
068  ENTER^
069  X<> IND 09
070  -
071  ABS
072  X<>Y
073  ENTER^
074  X<> IND 08
075  -
076  ABS
077  +
078  ST+ 13
079  ISG 10
080  ISG 08
081  ISG 09
082  GTO 02
083  RCL 13
084  VIEW X
085   E-9                     or another "small" number
086  X<Y?
087  GTO 01
088  RCL 06
089  ST- 11
090  ST- 12
091  CLST
092  LBL 04
093  X<>Y
094  RCL IND 10
095  RCL IND 11
096  *
097  +
098  X<>Y
099  RCL IND 10
100  RCL IND 12
101  *
102  +
103  ISG 10
104  ISG 11
105  ISG 12
106  GTO 04
107  ST+ 03
108  X<>Y
109  ST+ 02
110  RCL 04
111  ST+ 01
112  DSE 07
113  GTO 01           a three-byte GTO
114  RCL 03
115  RCL 02
116  RCL 01
117  CLD
118  END

( 177 bytes / SIZE n2+4n+14 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

-If, for instance, you are using the 6-stage 12-order method described above, the same 48 coefficients are to be stored into registers R26 thru R73

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

-Initialize registers R26 thru R73.

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03   and if we choose  h = 0.5
0.5  STO 04
2    STO 05    ( the number of steps )
6    STO 06    ( the number of stages )

XEQ "IRKB"   yields          x  =  1                           execution time = 10mn43s
RDN                         y(1) =  0.3678794411      ( error = -10-10 )
RDN                      z(0.5) = -0.7357588823      ( error = 0 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | k2n - k'2n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

2°)  Implicit Runge-Kutta methods, a special case:

A) First order differential equations:  dy/dx = f (x;y)

-Implicit Runge-Kutta methods are usually more "stable" than any explicit method, but n-stage (2n-1) or (2n-2) order methods are even more stable
than n-stage 2n-order methods. ( For a detailed description of stability properties of Runge-Kutta methods, cf Butcher's work )
-That's precisely what happens with the Radau IIA methods ( n-stage (2n-1)-order methods ) and Lobatto IIIC methods ( n-stage (2n-2)-order methods )
which also have another feature in common:  ci = an,i  ( i = 1 , 2 , .... , n ) , thus, less data registers are needed.
-The 2 following programs take advantage of this fact but the 2 programs presented above would work too.

Note:   The last programs listed in "Runge-Kutta methods for the HP-41" use a 5-stage 8-order Lobatto IIIC method.

Data Registers:     Registers R00 to R05   and    Rn+12 to Rn2+2n+11     are to be initialized before executing "I2RK"

•   R00:  f name   •   R01 = x0                                      R06 to R11:  temp             •  Rn+12 to Rn2+2n+11 =  the ( n2 + n ) coefficients ai,j ; bi
•   R02 = y0                                      R12 = k1                                                                     which are to be stored as shown below:
•   R03 = h = stepsize                        R13 = k2
•   R04 = N = number of steps          ............
•   R05 = n  = number of stages        Rn+11 = kn

•  Rn+12   = a1,1        •  Rn+13   = a1,2   ...........................    •  R2n+11 = a1,n       •  R2n+12 = b1
•  R2n+13 = a2,1        •  R2n+14 = a2,2   ...........................    •  R3n+12 = a2,n       •  R3n+13 = b2

..................................................................................................................................................

•  Rn2+n+11 = an,1     •  Rn2+n+12 = an,2   ........................   •  Rn2+2n+10 = an,n  •  Rn2+2n+11 = bn

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

01  LBL "I2RK"
02  RCL 04
03  STO 06
04  RCL 05
05  11
06  +
07  .1
08  %
09  12
10  +
11  STO 07
12  CLRGX               if you don't have an HP-41CX , replace line 12 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00
13  LBL 01
14  RCL 05
15  STO 08
16  RCL 07
17  FRC
18  12
19  ST+ 08
20  +
21  STO 07
22  CLX
23  ST0 10
24  LBL 02
25  RCL 07
26  FRC
27  12
28  +
29  STO 09
30  CLX
31  LBL 03
32  RCL IND 08
33  RCL IND 09
34  *
35  +
36  ISG 08
37  CLX
38  ISG 09
39  GTO 03
40  RCL 02
41  +
42  STO 11
43  RCL 03
44  RCL IND 08
45  *
46  RCL 01
47  +
48  XEQ IND 00
49  RCL 03
50  *
51  ENTER^
52  X<> IND 07
53  -
54  ABS
55  ST+ 10
56  ISG 08
57  CLX
58  ISG 07
59  GTO 02
60  RCL 10
61  VIEW X
62   E-9                     or another "small" number
63  X<Y?
64  GTO 01
65  RCL 03
66  ST+ 01
67  RCL 11
68  STO 02
69  DSE 06
70  GTO 01
71  RCL 01
72  CLD
73  END

( 109 bytes / SIZE n2+2n+12 )

 STACK INPUTS OUTPUTS Y / y(x0+N.h) X / x0+N.h

-For instance, here are the 56 coefficients of a 7-stage 13-order method ( a Radau IIA method )
-These coefficients are given with 20 decimals in case you want to use them with a more accurate calculator:

R19 = a1,1 =  0.03754 62649 93921 33133          R27 = a2,1 =  0.08014 75965 15618 96780          R35 = a3,1 =  0.07206 38469 41881 90211
R20 = a1,2 = -0.01403 93345 56460 40154          R28 = a2,2 =  0.08106 20639 85891 53668          R36 = a3,2 =  0.17106 83549 83886 61942
R21 = a1,3 =  0.01035 27896 00742 30094          R29 = a2,3 = -0.02123 79921 20711 03494          R37 = a3,3 =  0.10961 45640 40072 10923
R22 = a1,4 = -0.00815 83225 40275 01191          R30 = a2,4 =  0.01400 02912 38817 11898          R38 = a3,4 = -0.02461 98717 28984 05386
R23 = a1,5 =  0.00638 84138 79534 68494          R31 = a2,5 = -0.01023 41857 30090 16383          R39 = a3,5 =  0.01476 03770 43950 81707
R24 = a1,6 = -0.00460 23267 79148 65550          R32 = a2,6 =  0.00715 34651 51364 59050          R40 = a3,6 = -0.00957 52593 96791 40056
R25 = a1,7 =  0.00182 89425 61470 64370          R33 = a2,7 = -0.00281 26393 72406 72334          R41 = a3,7 =  0.00367 26783 97138 30567
R26 =  b1  =  0.02931 64271 59784 89197          R34 =  b2  =  0.14807 85996 68484 29185          R42 =   b3 =   0.33698 46902 81154 29910

R43 = a4,1 =  0.07570 51258 19824 42042          R51 = a5,1 =  0.07391 23421 63191 84654          R59 = a6,1 =  0.07470 55620 59796 23017
R44 = a4,2 =  0.15409 01551 42171 14465          R52 = a5,2 =  0.16135 56076 15942 43219          R60 = a6,2 =  0.15830 72238 72468 70066
R45 = a4,3 =  0.22710 77366 73202 38641          R53 = a5,3 =  0.20686 72415 52104 19782          R61 = a6,3 =  0.21415 34232 67200 03111
R46 = a4,4 =  0.11747 81870 37024 78199          R54 = a5,4 =  0.23700 71153 42694 23476          R62 = a6,4 =  0.21987 78470 31860 03999
R47 = a4,5 = -0.02381 08271 53044 17358         R55 = a5,5 =  0.10308 67935 33813 44662           R63 = a6,5 =  0.19875 21216 80635 26980
R48 = a4,6 =  0.01270 99855 33661 20563          R56 = a5,6 = -0.01885 41391 52580 44884          R64 = a6,6 =  0.06926 55016 05509 13323
R49 = a4,7 = -0.00460 88442 81289 63344         R57 = a5,7 =  0.00585 89009 74888 79182           R65 = a6,7 = -0.00811 60081 97728 29011
R50 =  b4  =  0.55867 15187 71550 13208          R58 =  b5  =  0.76923 38620 30054 50092           R66 =  b6  =  0.92694 56713 19741 11485

R67 = a7,1 =  0.07449 42355 56010 31793
R68 = a7,2 =  0.15910 21157 33650 74087
R69 = a7,3 =  0.21235 18895 02977 80420
R70 = a7,4 =  0.22355 49145 07283 23475                         (   and    ci = a7,i  )
R71 = a7,5 =  0.19047 49368 22115 57690
R72 = a7,6 =  0.11961 37446 12656 20289
R73 = a7,7 =  0.02040 81632 65306 12245 = 1/49
R74 =  b7  =  1

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

-Initialize registers R19 thru R74

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.5
0.5  STO 03
2    STO 04    ( the number of steps )
7    STO 05    ( we are using a 7-stage method )                      XEQ "I2RK"   yields          x = 1                         ( in 7mn13s )
X<>Y                  y(1) = 0.3678794410     ( error = -2 10-10 )

( With  h = N = 1 the error is only  -1.5 10-9 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | kn - k'n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

B) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:     Registers R00 to R06 and    R2n+14 to Rn2+4n+13     are to be initialized before executing "I2RKB"

•   R00:  f name   •   R01 = x0                                      R07 to R15:  temp             •  R2n+16 to Rn2+3n+15 =  the ( n2 + n ) coefficients ai,j ; bi
•   R02 = y0                                      R16 to R2n+15 = ki                                                      which are to be stored as shown below:
•   R03 = z0
•   R04 = h   = stepsize
•   R05 = N  = number of steps
•   R06 = n   = number of stages

•  R2n+16 = a1,1        •  R2n+17 = a1,2   ...........................    •  R3n+15 = a1,n         •  R3n+16 = b1
•  R3n+17 = a2,1        •  R3n+18 = a2,2   ...........................    •  R4n+16 = a2,n         •  R4n+17 = b2

.....................................................................................................................................................

•  Rn2+2n+15 = an,1   •  Rn2+2n+16 = an,2   ......................   •  Rn2+3n+14 = an,n   •  Rn2+3n+15 = bn

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

001  LBL "I2RKB"
002  RCL 05
003  STO 07
004  RCL 06
005  ST+ X
006  15
007  +
008  .1
009  %
010  16
011  +
012  STO 09
013  CLRGX               if you don't have an HP-41CX , replace line 13 by the 5 lines:      0    LBL 00    STO IND Y    ISG Y    GTO 00
014  LBL 01
015  RCL 09
016  FRC
017  16
018  STO 10
019  +
020  STO 08
021  RCL 06
022  ST+ 10
023  ST+ 10
024  +
025  STO 09
026  CLX
027  STO 13
028  LBL 02
029  RCL 09
030  FRC
031  16
032  +
033  STO 11
034  RCL 06
035  +
036  STO 12
037  CLST
038  LBL 03
039  X<>Y
040  RCL IND 10
041  RCL IND 11
042  *
043  +
044  X<>Y
045  RCL IND 10
046  RCL IND 12
047  *
048  +
049  ISG 10
050  CLX
051  ISG 11
052  ISG 12
053  GTO 03
054  RCL 03
055  +
056  STO 15
057  X<>Y
058  RCL 02
059  +
060  STO 14
061  RCL 04
062  RCL IND 10
063  *
064  RCL 01
065  +
066  XEQ IND 00
067  RCL 04
068  ST* Z
069  *
070  ENTER^
071  X<> IND 09
072  -
073  ABS
074  X<>Y
075  ENTER^
076  X<> IND 08
077  -
078  ABS
079  +
080  ST+ 13
081  ISG 10
082  CLX
083  ISG 08
084  ISG 09
085  GTO 02
086  RCL 13
087  VIEW X
088   E-9                     or another "small" number
089  X<Y?
090  GTO 01
091  RCL 04
092  ST+ 01
093  RCL 15
094  STO 03
095  RCL 14
096  STO 02
097  DSE 07
098  GTO 01               a three-byte GTO
099  RCL 01
100  CLD
101  END

( 146 bytes / SIZE  n2+3n+16 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

-If, for instance, you are using the 7-stage 13-order method described above, the same 56 coefficients are to be stored into registers R30 thru R85.

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

-Initialize registers R30 thru R85.

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03   and if we choose  h = 0.5
0.5  STO 04
2    STO 05    ( the number of steps )
7    STO 06    ( the number of stages )

XEQ "I2RKB"   yields        x  =  1                           execution time = 13mn54s
RDN                         y(1) =  0.3678794411      ( error = -10-10 )
RDN                      z(0.5) = -0.7357588822      ( error =  10-10 )

-The HP-41 displays  | k1 - k'1 | + ............ +  | k2n - k'2n |  where   ki and k'i  are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize h and start over!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.

Reference:     J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5