hp41programs

Runge-Kutta10

Runge-Kutta Formula of order 10 for the HP-41


Overview
 

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

-These routines apply a 17-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 = 17 so  169 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  )

    k17 = h.f(x+b17.h;y+a17,1.k1+a17,2.k2+ .......... +a17,16.k16)

then,     y(x+h) = y(x) + c1.k1+ ................ + c17.k17
 

0°)  169 Constants
 

-All the programs listed in paragraphs 1-2-3 use 169 constants.
-"C169" stores these numbers into registers R65 thru R233
-Registers R07 to R64 are also used.

-Allways execute "C169" before executing the other programs... unless you store your own 169 constants into R65 to R233 !
 

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

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

                                         •  R200= a17,1  ........................................    •  R215 = a17,16       •  R216 = b17

                                         •  R217 = c1 .........................................................................       •  R233 = c17
 
 
 
 

 01 LBL "C169"
 02 8.085
 03 CLRGX
 04 .8950802958
 05 STO 07
 06 .1979668312
 07 STO 09
 08 .0729547847
 09 CHS
 10 STO 10
 11 .8512362396
 12 CHS
 13 STO 12
 14 .3983201123
 15 STO 13
 16 3.639372632
 17 STO 14
 18 1.54822877
 19 STO 15
 20 2.122217147
 21 CHS
 22 STO 16
 23 1.583503985
 24 CHS
 25 STO 17
 26 1.715616082
 27 CHS
 28 STO 18
 29 .0244036406
 30 CHS
 31 STO 19
 32 .3090367612
 33 STO 20
 34 .9151765614
 35 CHS
 36 STO 21
 37 1.454534402
 38 STO 22
 39 .7773336436		  40 STO 34
 41 CHS
 42 STO 25
 43 .0910895662
 44 STO 33
 45 CHS
 46 STO 27
 47 .5393578408
 48 STO 35
 49 .1
 50 STO 36
 51 STO 51
 52 STO 61
 53 .1571786658
 54 STO 50
 55 CHS
 56 STO 38
 57 .1817813007
 58 STO 52
 59 .675
 60 STO 53
 61 CHS
 62 STO 67
 63 .3427581598
 64 STO 54
 65 CHS
 66 STO 66
 67 .2591112146
 68 STO 56
 69 CHS
 70 STO 65
 71 .3582789667
 72 CHS
 73 STO 57
 74 1.045948959
 75 CHS
 76 STO 58
 77 .9303278454
 78 STO 59		  79 1.779509594
 80 STO 60
 81 .2825475695
 82 CHS
 83 STO 62
 84 .1593273501
 85 CHS
 86 STO 63
 87 .1455158946
 88 CHS
 89 STO 64
 90 1
 91 STO 68
 92 30
 93 1/X
 94 STO 69
 95 STO 71
 96 STO 85
 97 CHS
 98 STO 83
 99 40
100 1/X
101 STO 70
102 CHS
103 STO 84
104 .05
105 STO 73
106 CHS
107 STO 82
108 .04
109 STO 75
110 CHS
111 STO 81
112 .1892374781
113 STO 77
114 STO 80
115 .2774291885
116 STO 78
117 STO 79		 118 7.155079
119 REGMOVE
120 9.096
121 CLRGX
122 .1
123 STO 07
124 STO 08
125 STO 83
126 .9151765614
127 CHS
128 STO 09
129 1.454534402
130 STO 10
131 .5393578408
132 STO 11
133 .2022591903
134 ENTER
135 STO 12
136 3
137 *
138 STO 14
139 +
140 STO 15
141 .5
142 -
143 STO 20
144 .1840247147
145 STO 16
146 .1979668312
147 STO 18
148 .0729547847
149 CHS
150 STO 19
151 .087900734
152 STO 21
153 .4104597025
154 STO 24
155 .4827137537
156 STO 25		 157 .9810741902
158 STO 26
159 .0859700505
160 STO 27
161 .330885963
162 STO 30
163 .4896629573
164 STO 31
165 .0731856375
166 CHS
167 STO 32
168 1.2
169 1/X
170 STO 33
171 STO 96
172 .1209304491
173 STO 34
174 .2601246758
175 STO 38
176 .0325402622
177 STO 39
178 .0595780212
179 CHS
180 STO 40
181 .3540173659
182 STO 41
183 .1108543796
184 STO 42
185 .0605761488
186 CHS
187 STO 47
188 .3217637056
189 STO 48
190 .5104857256
191 STO 49
192 .882527662
193 STO 50
194 .1120544148
195 STO 51		 196 .1449427759
197 CHS
198 STO 56
199 .3332697191
200 CHS
201 STO 57
202 .4992692296
203 STO 58
204 .5095046089
205 STO 59
206 .6426157582
207 STO 60
208 .113976784
209 STO 61
210 .0768813364
211 CHS
212 STO 66
213 .2395273603
214 STO 67
215 .3977746624
216 STO 68
217 .0107558957
218 STO 69
219 .3277691242
220 CHS
221 STO 70
222 .3573842418
223 STO 71
224 .0798314528
225 STO 72
226 .0520329687
227 CHS
228 STO 77
229 .0576954146
230 CHS
231 STO 78
232 .1947819157
233 STO 79		 234 .1453849232
235 STO 80
236 .078294271
237 CHS
238 STO 81
239 .1145032994
240 CHS
241 STO 82
242 .117472338
243 STO 83
244 .9851156101
245 STO 84
246 .330885963
247 STO 87
248 .4896629573
249 STO 88
250 1.378964866
251 CHS
252 STO 89
253 .861164195
254 CHS
255 STO 90
256 5.784288136
257 STO 91
258 3.28807762
259 STO 92
260 2.386339051
261 CHS
262 STO 93
263 3.254793425
264 CHS
265 STO 94
266 2.163435417
267 CHS
268 STO 95
269 7.06509
270 REGMOVE
271 END

 
     ( 1271 bytes / SIZE 234 )
 
 

      STACK        INPUT      OUTPUT
           X             /             /

 
-Execution time = 35 seconds.
 
 

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 R233 are to be initialized before executing "RK10" )

                                     •  R01 = x0      •  R03 = h  = stepsize                 R05 & R10: temp            •  R65 to •  R233:  the 169 coefficients stored by "C169"
                                     •  R02 = y0      •  R04 = N = number of steps    R11 to R27:  k1 to k17

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.027
 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 17
 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

 
     ( 117 bytes / SIZE 234 )
 
 

      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 "C169"

 "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 = 11m12s---
                              X<>Y                       y(1) = 0.3678794412

-The exact result is 1/e =  0.3678794412... so there is no error at all !
-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 R233 are to be initialized before executing "2RK10" )

                                     •  R01 = x0      •  R04 = h  = stepsize                 R06 & R10: temp            •  R65 to •  R233:  the 169 coefficients stored by "C169"
                                     •  R02 = y0      •  R05 = N = number of steps    R11 to R27:  k1 to k17
                                     •  R03 = z0                                                        R28 to R44:  l1 to l17
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 45
 06 15.999
 07 STO 46
 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 16
 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 45
 51 +
 52 STO 08
 53 17.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 46
 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 234 )
 
 

      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 "C169"

 "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 = 19m00s---
                                  RDN                      y(1) =  0.3678794412
                                  RDN                      z(1) = -0.7357588823

-The exact results are  y(1) = 1/e =  0.3678794412... & z(1) = -2/e = -0.7357588823...  again no error at all !
-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 R233 are to be initialized before executing "3RK10" )

                                     •  R01 = x0      •  R05 = h  = stepsize                 R07 & R10: temp            •  R65 to •  R233:  the 169 coefficients stored by "C169"
                                     •  R02 = y0      •  R06 = N = number of steps    R11 to R27:  k1 to k17         R62 to R64: temp
                                     •  R03 = z0                                                        R28 to R44:  l1 to l17
                                     •  R04 = t0                                                        R45 to R61:  m1 to m17
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 62
 04 GTO 01
 05 LBL 00
 06 RCL IND 07
 07 RCL IND 10
 08 *
 09 ST+ 64
 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 64
 31 RCL 02
 32 +
 33 RTN
 34 LBL 01
 35 16
 36 STO 63
 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 45
 51 RDN
 52 STO 28
 53 X<>Y
 54 STO 11
 55 LBL 02
 56 RCL 07
 57 FRC
 58 11.001
 59 +
 60 STO 07
 61 17.9
 62 +
 63 STO 08        
 64 LASTX
 65 INT
 66 +
 67 STO 09 
 68 CLST
 69 STO 64		  70 XEQ 00
 71 STO 64
 72 CLX
 73 RCL 05
 74 RCL IND 10
 75 *
 76 RCL 01
 77 +
 78 RCL 64
 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 63		  92 GTO 02
 93 15.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 62
108 GTO 01
109 RCL 04       
110 RCL 03 
111 RCL 02
112 RCL 01
113 END

 
     ( 189 bytes / SIZE 234 )
 
 

      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 "C169"

 "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 = 27m55s---
                                                     RDN     y(1) =  0.258207907
                                                     RDN     z(1) =  1.157623979
                                                     RDN     t(1) =  0.842178312

-The exact values are:

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

-Here, we get a ( small ) error.
-Press R/S to continue with the same h and N or modify   R05-R06

-Synthetic registers M , N , O , P , Q  and register R64 may also be used in the subroutine to calculate f , g , h
-If it's not enough, store the constants in - for example - R71 to R239 and replace line 39 ( 65.9 ) by 71.9
 

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

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

                                      •  R156+20n = a17,1  ..............................    •  R171+20n = a17,16       •  R172+20n = b17

                                      •  R173+20n = c1 ........................................................................       •  R189+20n = c17
 

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.018
 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 190+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

 
We store the name of the subroutine in R00    "T"  ASTO 00
the numbers of equations in R01                      3    STO 01

-Then we have to store the 169 constants into R21+20n to R189+20n
-Execute for instance the following routine:
 
 

 01 LBL "INIT"
 02 XEQ "C169"
 03 RCL 01
 04 20
 05 *
 06 21.169
 07 +
 08  E3
 09 /
 10 65
 11 +
 12 REGMOVE
 13 END

 
-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 = 5s---            with V41 in turbo mode
                                               RDN     y(1) =  0.258207907
                                               RDN     z(1) =  1.157623977
                                               RDN     t(1) =  0.842178313

Notes:

-With a real HP41, execution time is about 50 minutes, but with free42, the results are returned almost instantaneously !
-Since SIZE = 190+20n, "NRK10" can solve at most a system of n = 6 differential equations with SIZE 310.
-But with free42, n may reach 40 with SIZE 990 ( 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]  T. Feagin - "A Tenth-Order Runge-Kutta Method with Error Estimate"

-These programs employ the coefficients given in reference [2] without taking into account the error estimate.