hp41programs

Runge-Kutta10B

Runge-Kutta Formula of order 10 (2) for the HP-41


Overview
 

 0°)  151 Constants
 1°)  1 Differential Equation
 2°)  System of 2 Differential Equations
 3°)  System of 3 Differential Equations
 4°)  System of n Differential Equations
 

-Thanks to Valentin Albillo, I've found on the web a 16-stage 10th-order Runge-Kutta formula to solve ODEs

-All s-stage explicit Runge-Kutta formulae without error estimate are determined by ( s2 + 3.s - 2 ) / 2  coefficients ai,j ; bi ; ci
-Here, s = 16 so there are 151 constants.

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

    k16 = h.f(x+b16.h;y+a16,1.k1+a16,2.k2+ .......... +a16,15.k15)

then,     y(x+h) = y(x) + c1.k1+ ................ + c16.k16
 

0°)  151 Constants
 

-All the programs listed in paragraphs 1-2-3 use 151 constants.
-"C151" stores these numbers into registers R65 thru R215
-Registers R01 to R99 are also used.

-Allways execute "C151" before executing the other programs... unless you store your own 151 constants into R65 to R215
 

                                         •  R65 = a2,1                                                                                •  R66 = b2
                                         •  R67 = a3,1        •  R68 = a3,2                                                    •  R69 = b3

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

                                         •  R184= a16,1  ........................................    •  R198 = a16,15       •  R199 = b16

                                         •  R200 = c1 .........................................................................       •  R215 = c16
 
 
 
 

 01 LBL "C151"
 02 .0688809661
 03 STO 01
 04 STO 02
 05 STO 09
 06 -.8381052035
 07 STO 03
 08 1.253690049
 09 STO 04
 10 .4155848452
 11 STO 05
 12 -.0049094868
 13 STO 06
 14 .0823217303
 15 STO 07
 16 -.0085312774
 17 STO 08
 18 1.048931954
 19 STO 10
 20 -.7538281773
 21 STO 11
 22 .8052281597
 23 STO 12
 24 -.3844632207
 25 STO 13
 26 .7158687154
 27 STO 14
 28 -.2399238343
 29 STO 15
 30 -.0636426116
 31 STO 16
 32 .1996754314
 33 STO 17
 34 .6103537951
 35 STO 18
 36 .3785994723
 37 STO 19
 38 .8850622527
 39 STO 20
 40 .0177883395
 41 STO 21
 42 -.0110502191
 43 STO 22
 44 -.0043934286
 45 STO 23
 46 .1059752729
 47 STO 24
 48 .0040508907
 49 STO 25
 50 -.0016792971
 51 STO 26
 52 .1106915584
 53 STO 27
 54 .2356604623
 55 STO 28
 56 .0889336269
 57 STO 29
 58 .0413883471
 59 STO 30
 60 -.852903036
 61 STO 31
 62 -.0230808755
 63 STO 32
 64 .0096859219
 65 STO 33
 66 .8014497793
 67 STO 34
 68 .301134226
 69 STO 35
 70 .0904869376
 71 STO 36
 72 .032044272
 73 STO 37
 74 .1243376822
 75 STO 38
 76 -.3073152176
 77 STO 39
 78 .1644784368
 79 STO 40
 80 -.0410068117
 81 STO 41
 82 .4066198993
 83 STO 42
 84 .1866998712
 85 STO 43
 86 .6563450699
 87 STO 44
 88 -.1281584914
 89 STO 45
 90 -.0949429224
 91 STO 46
 92 -.1445034451
 93 STO 47
 94 .921349407
 95 STO 48
 96 -.1326530105
 97 STO 49
 98 .0166148587
 99 STO 50
100 -.6446177324
101 STO 51
102 .4714846368
103 STO 52
104 .1510115445
105 STO 53
106 .4155848452
107 STO 54
108 STO 65
109 -.3539426289
110 STO 55
111 -.1299925006
112 STO 56
113 .0729671747
114 STO 57
115 1.234092882
116 STO 58
117 .2432547931
118 STO 59
119 -.0488641534
120 STO 60
121 -.5144237141
122 STO 61
123 .0708850045
124 STO 62
125 -.2055533399
126 STO 63
127 .0471613277
128 STO 64
129 -1.288835964
130 STO 66
131 -.2802086976
132 STO 67
133 2.955025311
134 STO 68
135 2.687294528
136 STO 69
137 4.94062343
138 STO 70
139 -1.10669748
140 STO 71
141 -1.014991843
142 STO 72
143 2.537509826
144 STO 73
145 -3.050854355
146 STO 74
147 -4.48248369
148 STO 75
149 -1.212578432
150 STO 76
151 .6838026329
152 STO 77
153 -.2644034429
154 STO 78
155 -.0668111881
156 STO 79
157 .2298317881
158 STO 80
159 .6407414965
160 STO 81
161 .4294291108
162 STO 82
163 -.0096707888
164 STO 83
165 .0161366667
166 STO 84
167 -.0359631432
168 STO 85
169 -.0490267009
170 STO 86
171 .0138853767
172 STO 87
173 -.0232947511
174 STO 88
175 .0041107147
176 STO 89
177 .8849651386
178 STO 90
179 1.289959279
180 STO 91
181 .135032652
182 STO 92
183 -1.575337353
184 STO 93
185 -1.136141441
186 STO 94
187 -2.606299763
188 STO 95
189 -.1422551486
190 STO 96
191 -1.300502794
192 STO 97
193 2.707971631
194 STO 98
195 2.435542969
196 STO 99
197 1.065099
198 REGMOVE
199 -.2645818252
200 STO 01
201 .2843174728
202 STO 02
203 -.015673365
204 STO 03
205 .6035525316
206 STO 04
207 .4155848452
208 STO 05
209 STO 20
210 -.8217014524
211 STO 06
212 1.255912216
213 STO 07
214 -.0161528412
215 STO 08
216 -.0229703899
217 STO 09
218 -.0279603438
219 STO 10
220 -.7509429873
221 STO 11
222 -.0083102909
223 STO 12
224 .0278552386
225 STO 13
226 .0402142154
227 STO 14
228 1.617136464
229 STO 15
230 -1.396283777
231 STO 16
232 -.0142443916
233 STO 17
234 .7570905434
235 STO 18
236 -.2240573576
237 STO 19
238 .2772640185
239 STO 21
240 .0945381837
241 STO 22
242 .841727543
243 STO 23
244 -.9066525983
245 STO 24
246 -.0933485803
247 STO 25
248 4.088877589
249 STO 26
250 .795399865
251 STO 27
252 -.0485917515
253 STO 28
254 .1481985146
255 STO 29
256 -1.770211481
257 STO 30
258 1.838210845
259 STO 31
260 .0490934446
261 STO 32
262 -3.803788231
263 STO 33
264 .3315723987
265 STO 34
266 -.8422897608
267 STO 35
268 1
269 STO 36
270 .0318192746
271 STO 37
272 CLX
273 STO 38
274 STO 40
275 STO 41
276 .0468136929
277 STO 39
278 1.375535362
279 STO 42
280 .1750656714
281 STO 43
282 .1492479865
283 STO 44
284 .2718099231
285 STO 45
286 -.1707794051
287 STO 46
288 .3003551977
289 STO 47
290 -.010142544
291 STO 48
292 -1.191778158
293 STO 49
294 .0363913935
295 STO 50
296 -.0468331609
297 STO 51
298 .0324947663
299 STO 52
300 1.164052
301 REGMOVE
302 END

 
     ( 1929 bytes / SIZE 216 )
 
 

      STACK        INPUT      OUTPUT
           X             /             /

 
-Execution time = 55 seconds.
 

Note:

-There are of course less constants than with a 17-stage Runge-Kutta method,
  but there are much more non-zero coefficients: so "C151" is longer & slower than "C169" !
 

1°)  1 Differential Equation
 

-We want to solve  y' = f(x,y)  with the initial conditions    y(x0) = y0
 
 

Data Registers:          •  R00:  f name                                     ( Registers R00 to R04 & R65 to R215 are to be initialized before executing "RK10" )

                                     •  R01 = x0      •  R03 = h  = stepsize                 R05 & R10: temp            •  R65 to •  R215:  the 151 coefficients stored by "C151"
                                     •  R02 = y0      •  R04 = N = number of steps    R11 to R26:  k1 to k16

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

 01 LBL "RK10" 
 02 RCL 04
 03 STO 06
 04 12.026
 05 STO 05
 06 11.001
 07 STO 10
 08 LBL 01
 09 .01
 10 STO 07
 11 65.9
 12 STO 08
 13 RCL 05
 14 STO 09
 15 RCL 02
 16 RCL 01
 17 XEQ IND 00
 18 RCL 03
 19 *
 20 STO 11
 21 LBL 02
 22 RCL 07
 23 FRC
 24 RCL 10
 25 +
 26 STO 07
 27 CLX
 28 LBL 03
 29 RCL IND 07
 30 RCL IND 08
 31 *
 32 +
 33 ISG 08
 34 ISG 07
 35 GTO 03
 36 RCL 02
 37 +
 38 RCL 03
 39 RCL IND 08
 40 *
 41 RCL 01
 42 +
 43 XEQ IND 00
 44 RCL 03
 45 *
 46 STO IND 09
 47 ISG 08
 48 ISG 09
 49 GTO 02
 50 16
 51 ST- 09
 52 CLX
 53 LBL 04
 54 RCL IND 08
 55 RCL IND 09
 56 *
 57 +
 58 ISG 08
 59 ISG 09
 60 GTO 04 
 61 ST+ 02
 62 RCL 03
 63 ST+ 01
 64 DSE 06
 65 GTO 01
 66 RCL 02        
 67 RCL 01
 68 END

 
     ( 115 bytes / SIZE 216 )
 
 

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

 
Example:    y' = -2 x.y   y(0) = 1  y'(0) = 0      Evaluate  y(1)
 
 

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

 
   XEQ "C151"

 "T"  ASTO 00

  0     STO 01
  1     STO 02    and if we choose  h = 0.1
 0.1   STO 03
 10    STO 04    XEQ "RK10"   yields          x = 1                                       ---Execution time = 10m01s---
                              X<>Y                       y(1) = 0.3678794410

-The exact result is 1/e =  0.3678794412...
-Press R/S to continue with the same h and N or modify R03-R04
 

2°)  System of 2 Differential Equations
 

-"2RK10" solves the system

    y' = f(x,y,z)      with the initial values     y(x0) = y0
    z' = g(x,y,z)                                         z(x0) = z0
 
 

Data Registers:          •  R00:  f name                                   ( Registers R00 to R05 & R65 to R215 are to be initialized before executing "2RK10" )

                                     •  R01 = x0      •  R04 = h  = stepsize                 R06 & R10: temp            •  R65 to •  R215:  the 151 coefficients stored by "C151"
                                     •  R02 = y0      •  R05 = N = number of steps    R11 to R26:  k1 to k16
                                     •  R03 = z0                                                        R27 to R42:  l1 to l16
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.
 
 
 
 

 01 LBL "2RK10"
 02 RCL 05
 03 STO 07
 04 11.001
 05 STO 43
 06 14.999
 07 STO 44
 08 GTO 01
 09 LBL 00
 10 X<>Y
 11 RCL IND 08
 12 RCL IND 10
 13 *
 14 +
 15 X<>Y
 16 RCL IND 09
 17 RCL IND 10
 18 *
 19 +
 20 ISG 10
 21 ISG 09
 22 ISG 08
 23 GTO 00
 24 RCL 03 
 25 +
 26 X<>Y
 27 RCL 02
 28 +
 29 RTN
 30 LBL 01 
 31 .01
 32 STO 08        
 33 15
 34 STO 06
 35 65.9
 36 STO 10
 37 RCL 03
 38 RCL 02
 39 RCL 01
 40 XEQ IND 00
 41 RCL 04
 42 ST* Z
 43 *
 44 STO 28
 45 X<>Y
 46 STO 11
 47 LBL 02
 48 RCL 08
 49 FRC
 50 RCL 43
 51 +
 52 STO 08
 53 16.9
 54 +
 55 STO 09
 56 CLST
 57 XEQ 00
 58 RCL 04
 59 RCL IND 10
 60 *
 61 RCL 01
 62 +
 63 XEQ IND 00
 64 RCL 04
 65 ST* Z
 66 *
 67 STO IND 09
 68 X<>Y
 69 STO IND 08
 70 ISG 10
 71 DSE 06
 72 GTO 02
 73 RCL 44
 74 ST- 08
 75 ST- 09
 76 CLST
 77 XEQ 00
 78 STO 02
 79 X<>Y
 80 STO 03
 81 RCL 04
 82 ST+ 01
 83 DSE 07
 84 GTO 01
 85 RCL 03 
 86 RCL 02        
 87 RCL 01
 88 END

 
      ( 150 bytes / SIZE 216 )
 
 

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

 
Example:     y(0) = 1 ; y'(0) = 0 ;  y'' + 2.x.y' + 2.y = 0   Evaluate y(1) and y'(1)

-This second order equation is equivalent to the system:             y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    where z = y'
 
 

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

 
  XEQ "C151"

 "U"  ASTO 00

  0      STO 01
  1      STO 02
  0      STO 03   and if we choose  h = 0.1
 0.1    STO 04
 10     STO 05    XEQ "2RK10"   yields       x =  1                                               ---Execution time = 17m01s---
                                  RDN                      y(1) =  0.3678794411
                                  RDN                      z(1) = -0.7357588821

-The exact results are  y(1) = 1/e =  0.3678794412... & z(1) = -2/e = -0.7357588823...
-Press R/S to continue with the same h and N or modify  R04-R05
 

3°)  System of 3 Differential Equations
 

-"3RK10" solves the system of 3 ODEs:

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

Data Registers:          •  R00:  f name                                   ( Registers R00 to R06 & R65 to R215 are to be initialized before executing "3RK10" )

                                     •  R01 = x0      •  R05 = h  = stepsize                 R07 & R10: temp            •  R65 to •  R215:  the 151 coefficients stored by "C151"
                                     •  R02 = y0      •  R06 = N = number of steps    R11 to R26:  k1 to k16         R59 to R61: temp
                                     •  R03 = z0                                                        R27 to R42:  l1 to l16
                                     •  R04 = t0                                                        R43 to R58:  m1 to m16
Flags: /
Subroutine:  A program that takes  x , y , z , t  in X , Y , Z , T and returns  f(x,y,z,t) , g(x,y,z,t) , h(x,y,z,t)  in registers  Z , Y , X  respectively.
 

-Line 108 is a three-byte GTO 01
 
 

 01 LBL "3RK10"
 02 RCL 06
 03 STO 59
 04 GTO 01
 05 LBL 00
 06 RCL IND 07
 07 RCL IND 10
 08 *
 09 ST+ 61
 10 X<> Z
 11 RCL IND 08
 12 RCL IND 10
 13 *
 14 +
 15 X<>Y
 16 RCL IND 09
 17 RCL IND 10
 18 *
 19 +
 20 ISG 10
 21 ISG 09
 22 ISG 08
 23 ISG 07
 24 GTO 00
 25 RCL 04
 26 +
 27 X<>Y
 28 RCL 03
 29 +
 30 RCL 61
 31 RCL 02
 32 +
 33 RTN
 34 LBL 01
 35 15
 36 STO 60
 37 .01
 38 STO 07
 39 65.9
 40 STO 10
 41 RCL 04
 42 RCL 03
 43 RCL 02
 44 RCL 01
 45 XEQ IND 00
 46 RCL 05
 47 ST* T
 48 ST* Z
 49 *
 50 STO 43
 51 RDN
 52 STO 27
 53 X<>Y
 54 STO 11
 55 LBL 02
 56 RCL 07
 57 FRC
 58 11.001
 59 +
 60 STO 07
 61 16.9
 62 +
 63 STO 08 
 64 LASTX
 65 INT
 66 +
 67 STO 09        
 68 CLST
 69 STO 64
 70 XEQ 00
 71 STO 61
 72 CLX
 73 RCL 05
 74 RCL IND 10
 75 *
 76 RCL 01
 77 +
 78 RCL 61
 79 X<>Y
 80 XEQ IND 00
 81 RCL 05
 82 ST* T
 83 ST* Z
 84 *
 85 STO IND 09
 86 RDN
 87 STO IND 08
 88 X<>Y
 89 STO IND 07
 90 ISG 10
 91 DSE 60
 92 GTO 02
 93 14.999
 94 ST- 07
 95 ST- 08
 96 ST- 09
 97 CLST
 98 STO 64
 99 XEQ 00
100 STO 02
101 RDN
102 STO 03
103 X<>Y
104 STO 04
105 RCL 05
106 ST+ 01
107 DSE 59
108 GTO 01
109 RCL 04 
110 RCL 03       
111 RCL 02
112 RCL 01
113 END

 
     ( 189 bytes / SIZE 216 )
 
 

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

 
Example:

     dy/dx = -y.z.t                    y(0) = 1         Evaluate  y(1) , z(1) , t(1)  with  h = 0.1   ( whence N = 10 )
     dz/dx = x.( y + z - t )         z(0) = 1
     dt/dx = x.y - z.t                  t(0) = 2
 
 

 01 LBL "V"
 02 R^
 03 CHS
 04 STO L
 05 R^
 06 ST* Y
 07 ST+ L
 08 CLX
 09 RCL Z
 10 R^
 11 ST+ L
 12 ST* Y
 13 X<> Z
 14 ST+ Y
 15 ST* Z
 16 X<> T
 17 LASTX
 18 *
 19 X<>Y
 20 RTN
 21 END

 
  XEQ "C151"

 "V"   ASTO 00

  0     STO 01
  1     STO 02  STO 03
  2     STO 04
 0.1   STO 05
 10    STO 06    XEQ "3RK10"   >>>>      x   =  1                                                       ---Execution time = 25m11s---
                                                     RDN     y(1) =  0.258207906
                                                     RDN     z(1) =  1.157623979
                                                     RDN     t(1) =  0.842178313

-The exact values are:

      y(1) = 0.258207906                      ( rounded to 9 decimals )
      z(1) = 1.157623981
      t(1) = 0.842178312

-Press R/S to continue with the same h and N or modify   R05-R06
-Synthetic registers M , N , O , P , Q  and registers R61-R62-R63-R64 may also be used in the subroutine to calculate f , g , h
 ( though register R61 is also used in "3RK10" )
 

Notes:

-These routines are slow with a real HP41
-But with V41 or 41CL, the results are obtained in a few seconds.
 

4°)  System of n Differential Equations
 

-Here, we want to solve the system:

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

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

Data Registers:           •  R00 = subroutine name

                     ( Registers R00 thru R03 , R20 thru R20+n and R21+20n thru R189+20n are to be initialized before executing "NRK10" )

                                      •  R01 = n  = number of equations = number of functions         R04 thru R10 & R19: temp    R11 to R18: unused
                                      •  R02 = h  = stepsize
                                      •  R03 = N = number of steps

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

                                      •  R20+n = yn(x0)
 

                                      •  R21+20n = a2,1                                                                                •  R22+20n = b2
                                      •  R23+20n = a3,1        •  R24+20n = a3,2                                            •  R25+20n = b3

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

                                      •  R140+20n = a16,1  ..............................    •  R154+20n = a16,15       •  R155+20n = b16

                                      •  R156+20n = c1 ........................................................................       •  R171+20n = c16
 

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 "NRK10"
 02 RCL 03
 03 STO 10
 04 RCL 01
 05 .1
 06 %
 07 +
 08 STO 05
 09 FRC
 10 21.02
 11 +
 12 STO 04
 13 RCL 05
 14 ST+ X
 15 +
 16 STO 06
 17 GTO 01
 18 LBL 02
 19 XEQ IND 00
 20 RCL 04
 21 RCL 05
 22 +
 23 RCL 07
 24 LBL 00
 25 RCL IND Y
 26 RCL 02
 27 *
 28 STO IND Y
 29 ISG Y
 30 RDN
 31 ISG Y
 32 GTO 00
 33 RTN
 34 LBL 03
 35 RCL IND Y
 36 STO IND Y
 37 ISG Y
 38 RDN
 39 ISG Y
 40 GTO 03
 41 RTN
 42 LBL 07
 43 RCL IND 07 
 44 RCL IND 08 
 45 *
 46 ST+ IND 04
 47 ISG 04
 48 TEXT0
 49 ISG 07
 50 GTO 07
 51 RCL 01
 52 ST- 04
 53 RCL 05
 54 FRC
 55 ST+ 07
 56 ISG 08
 57 GTO 07
 58 RTN
 59 LBL 01
 60 RCL 05
 61 RCL 06
 62 +
 63 STO 07
 64 XEQ 02
 65 RCL 04         
 66 RCL 06 
 67 XEQ 03
 68 RCL 01
 69 20
 70 *
 71 21
 72 +
 73 .1
 74 %
 75 +
 76 STO 08
 77 3.017
 78 STO 09
 79 RCL 20
 80 STO 19
 81 LBL 04
 82 XEQ 07
 83 RCL 02
 84 RCL IND 08 
 85 *
 86 ST+ 20
 87 XEQ 02
 88 RCL 06
 89 RCL 04
 90 XEQ 03
 91 RCL 19
 92 STO 20
 93 RCL 09
 94 INT
 95  E3
 96 /
 97 ISG X
 98 ST+ 08
 99 RCL 04
100 RCL 05
101 3
102 *
103 +
104 STO 07
105 ISG 09
106 GTO 04
107 RCL 02
108 ST+ 20
109 XEQ 07
110 DSE 10
111 GTO 01
112 RCL 23 
113 RCL 22        
114 RCL 21
115 RCL 20
116 END

 
     ( 191 bytes / SIZE 172+20.n )
 
 

      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:     Let's solve again 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 R20 ; R21 ; R22 ; R23 respectively, and we have to write a program to compute the 3 functions
 and store the results in R24 ; R25 ; R26   ( just after the "initial-value registers" ).

-For instance, the following subroutine:
 
 

 01  LBL "T"
 02  RCL 21
 03  RCL 22
 04  RCL 23
 05  *
 06  *
 07  CHS
 08  STO 24
 09  RCL 21
 10  RCL 22
 11  +
 12  RCL 23
 13  -
 14  RCL 20
 15  *
 16  STO 25
 17  RCL 20
 18  RCL 21
 19  *
 20  RCL 22
 21  RCL 23
 22  *
 23  -
 24  STO26
 25  RTN

 
-Then we have to store the 151 constants into R21+20n to R171+20n

  SIZE 232   XEQ "C151"  and since  n = 3

  65.081151   REGMOVE

-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 R20                      0   STO 20
     y1   in R21                      1   STO 21
     y2   in R22                      1   STO 22
     y3   in R23                      2   STO 23

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

    h is stored in R02                 0.1   STO 02
   N is stored in R03                 10    STO 03

             and  XEQ "NRK10"  >>>>        x   =  1                                                ---Execution time = 40m12s---
                                               RDN     y(1) =  0.258207906
                                               RDN     z(1) =  1.157623981
                                               RDN     t(1) =  0.842178312

Notes:

-With a real HP41, this program is very slow, but with free42, the results are returned almost instantaneously !
-Since SIZE = 172+20n, "NRK10" can solve at most a system of n = 7 differential equations with SIZE 312.
-But with free42, n may reach 41 with SIZE 992 ( beyond register 999, the control numbers would not work )
 
 

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  David K. Zhang - "Discovering New Runge-Kutta Methods Using Unstructured Numerical Search"