# hp41programs

Runge-Kutta-Nystrom10

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

Overview

0°)  116 Constants
1°)  1 Second Order Conservative Equation
2°)  System of 2 Second Order Conservative Equations
3°)  System of 3 Second Order Conservative Equations
4°)  System of n Second Order Conservative Equations

-These routines apply a 13-stage 10th-order Runge-Kutta-Nystrom formula to solve ODEs of the form:   y" = f(x,y)

-All s-stage explicit Runge-Kutta-Nystrom formulae without error estimate are determined by ( s2 + 5.s - 2 ) / 2 coefficients ai,j ; bi ; ci

-Here, s =13 , so 116 constants

k1 = h.f(x;y)
k2 = h.f(x+b2.h;y+b2.h.y'+h a2,1.k1)
k3 = h.f(x+b3.h;y+b3.h.y'+h(a3,1.k1+a3,2.k2) )
.................................................................................

k13 = h.f(x+b13.h;y+b13.h.y'+h (a13,1.k1+a13,2.k2+ .......... +a13,12.k12) )

then,     y(x+h) = y(x) + h y'(x) + h [ c1.k1+ ................ + c13.k13 ]
and      y'(x+h) = y'(x) + [ c'1.k1+ ................ + c'13.k13 ]

0°)  116 Constants

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

-Allways execute "C116" before executing the other programs... unless you store your own 116 constants into R65 to R180 !

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

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

•  R142= a13,1  ........................................    •  R153 = a13,12       •  R154 = b13

•  R155 = c1 .........................................................................       •  R167 = c13
•  R168 = c'1 ........................................................................       •  R180 = c'13

 01 LBL "C116"  02 7.068  03 CLRGX  04 .0507388591  05 CHS  06 STO 07  07 .852583158  08 CHS  09 STO 09  10 .2925628199  11 STO 10  12 .4263130454  13 CHS  14 STO 11  15 .3084512679  16 STO 12  17 .8206806306  18 STO 13  19 .2635320056  20 STO 14  21 .0380029595  22 CHS  23 STO 15  24 .05083695  25 STO 16  26 .8584004334  27 STO 17  28 .855063419  29 CHS  30 STO 18  31 4.417296783  32 CHS  33 STO 20 34 1.462041474  35 STO 21  36 1.580060367  37 STO 22  38 2.150047309  39 CHS  40 STO 23  41 5.219288295  42 STO 24  43 .701222246  44 CHS  45 STO 25  46 .4067428472  47 STO 26  48 .109950164  49 CHS  50 STO 27  51 .0254989511  52 STO 28  53 .9592205307  54 STO 29  55 3.612420577  56 STO 30  57 19.61437642  58 STO 32  59 6.554095453  60 CHS  61 STO 33  62 5.363477518  63 CHS  64 STO 34  65 8.954920063  66 STO 35 67 22.11199958  68 CHS  69 STO 36  70 3.066641833  71 STO 37  72 1.297438034  73 CHS  74 STO 38  75 .5982268483  76 STO 39  77 .0241070789  78 CHS  79 STO 40  80 .0045319136  81 STO 41  82 1  83 STO 42  84 .0114450454  85 STO 43  86 .1518222881  87 STO 48  88 .0938333833  89 STO 49  90 .131388714  91 STO 50  92 .042452794  93 STO 51  94 .0466143591  95 STO 52  96 .0197393408  97 STO 53  98 .0027040753  99 STO 54 100 .0114450454 101 STO 56 102 .1896345577 103 STO 61 104 .0990150484 105 STO 62 106 .2237772082 107 STO 63 108 .1046681151 109 STO 64 110 .1535817253 111 STO 65 112 .1394025506 113 STO 66 114 .0663097234 115 STO 67 116 .0121660259 117 STO 68 118 7.119062 119 REGMOVE 120 7.06 121 CLRGX 122 .0006255035 123 STO 07 124 .0353695786 125 STO 08 126 .0008340047 127 STO 09 128 .0016680095 129 STO 10 130 .0707391571 131 STO 11 132 .0143775922 133 STO 12 134 .0255203896 135 CHS 136 STO 13 137 .0299554191 138 STO 14 139 .1939722747 140 STO 15 141 .005660601 142 STO 16 143 .0224381472 144 STO 18 145 .009619895 146 STO 19 147 .2746584906 148 STO 20 149 .0042497841 150 STO 21 151 .0137373896 152 STO 23 153 .0021090824 154 STO 24 155 .0002169817 156 CHS 157 STO 25 158 .1993954583 159 STO 26 160 .0008446992 161 STO 27 162 .0007628617 163 STO 29 164 .004527489 165 CHS 166 STO 30 167 .0000947045 168 CHS 169 STO 31 170 .0043839568 171 STO 32 172 .0523320971 173 STO 33 174 .0402571887 175 STO 34 176 .2813255857 177 STO 36 178 .0990869733 179 CHS 180 STO 37 181 .0315586786 182 STO 38 183 .0762630093 184 STO 39 185 .2450911018 186 CHS 187 STO 40 188 .4128592672 189 STO 41 190 .5185220615 191 CHS 192 STO 42 193 4.022049675 194 CHS 195 STO 44 196 1.333954239 197 STO 45 198 .3628013762 199 CHS 200 STO 46 201 .446629326 202 CHS 203 STO 47 204 4.10933539 205 STO 48 206 .0833718601 207 STO 49 208 .5944056701 209 STO 50 210 .4552651504 211 STO 51 212 3.44102443 213 STO 53 214 1.142030463 215 CHS 216 STO 54 217 .3312250453 218 STO 55 219 .5251081451 220 STO 56 221 3.400720497 222 CHS 223 STO 57 224 .0149590255 225 STO 58 226 .017714834 227 STO 59 228 .6964849889 229 STO 60 230 7.065054 231 REGMOVE 232 END

( 1302 bytes / SIZE 181 )

 STACK INPUT OUTPUT X / /

---Execution time = 37s---

1°)  1 Differential Equation

-We want to solve  y" = f(x,y)

with the initial conditions    y(x0) = y0 ,  y'(x0) = y'0

Data Registers:          •  R00:  f name                                      ( Registers R00 to R05 & R65 to R180 are to be initialized before executing "RKN10" )

•  R01 = x0      •  R04 = h  = stepsize                 R06 & R09: temp            •  R65 to •  R180:  the 116 coefficients stored by "C116"
•  R02 = y0      •  R05 = N = number of steps    R11 to R23:  k1 to k13
•  R03 = y'0

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

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

( 138 bytes / SIZE 181 )

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

Example:     y(0) = 1  ,   y'(0) = 0  ,   y" = - y ( x2 + y2 ) 1/2

>>> Calculate  y(1) & y'(1)

 01  LBL "T"  02  X^2  03  X<>Y  04  STO Z  05  X^2  06  +  07  SQRT  08  *  09  CHS  10  RTN

"T"  ASTO 00

0  STO 01  STO 03  1  STO 02

-With h = 0.1  and  N = 10  ,  0.1  STO 04  10  STO 05

XEQ "RKS1"  >>>>      x   = 1                                            ---Execution time = 8m30s---
RDN   y(1) =  0.536630617
RDN   y'(1) = -0.860171927

Note:

-Simply press R/S to continue with the same h & N or modify registers R04 & R05 if need be.

2°)  System of 2 Differential Equations

-"2RKN10" 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:  f name                        ( Registers R00 to R07 & R65 to R180 are to be initialized before executing "2RKN10" )

•  R01 = x0      •  R04 = h  = stepsize                 •  R06 = y'0             R11 to R23:  k1 to k13     R08 to R10 & R37: temp
•  R02 = y0      •  R05 = N = number of steps    •  R07 = z'0             R24 to R36:  l1 to l13
•  R03 = z0                      •  R65 to •  R180:  the 116 coefficients stored by "C116"

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

( 180 bytes / SIZE 181 )

 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

 01  LBL "T"  02  RDN  03  STO Z  04  X<>Y  05  CHS  06  ST* Z  07  -  08  R^  09  *  10  RTN

"T"  ASTO 00

0     STO 01              1  STO 06
2     STO 02              1  STO 07
1     STO 03
0.1   STO 04
10    STO 05

XEQ "2RKN10"  >>>>     x   =  1                                                                                     ---Execution time = 13m58s---
RDN   y(1) = 1.531356645     and   R06 = y'(1) = -2.312840138
RDN   z(1) = 2.620254282              R07 = z'(1) =  2.941748401

Note:

-Simply press R/S to continue with the same h & N or modify registers R04 & R05 if need be.

3°)  System of 3 Differential Equations

-"3RKN10" solves the system:

y" = f(x,y,z,t)         y(x0) = y0     y'(x0) = y'0
z" = g(x,y,z,t)         z(x0) = z0     z'(x0) = z'0
t" = h(x,y,z,t)          t(x0) = t0      t'(x0) = t'0

Data Registers:          •  R00:  f name                        ( Registers R00 to R09 & R65 to R180 are to be initialized before executing "3RKN10" )

•  R01 = x0      •  R05 = h  = stepsize                 •  R07 = y'0             R16 to R28:  k1 to k13     R10 to R15 & R55: temp
•  R02 = y0      •  R06 = N = number of steps    •  R08 = z'0             R29 to R41:  l1 to l13
•  R03 = z0                                                        •  R09 = t'0             R42 to R54:  m1 to m13
•  R04 = t0                                     •  R65 to •  R180:  the 116 coefficients stored by "C116"

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

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

( 242 bytes / SIZE 181 )

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

Example:

d2y/dx2 = -y.z.t                    y(0) = 1        y'(0) = 1         Evaluate  y(1) , z(1) , t(1)  ,   y'(1) , z'(1) , t'(1)      with  h = 0.1  &  N = 10
d2z/dx2 = x.( y + z - t )         z(0) = 1        z'(0) = 1
d2u/dx2 = x.y - z.t                 t(0) = 2        t'(0) = 1

 01 LBL "T"  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

"T"   ASTO 00

0     STO 01
1     STO 02  STO 03  STO 07  STO 08  STO 09
2     STO 04
0.1   STO 05
10    STO 06    XEQ "3RKN10"  >>>>    x   = 1                                                                                          ---Execution time = 20m16s---

RDN   y(1) = 0.439524100           y'(1) = -2.101122880 = R07
RDN   z(1) = 2.070940654           z'(1) =   1.269596951 = R08
RDN   t(1) = 1.744524962            t'(1) = -1.704234756 = R09

Notes:

-The errors are about E-9 ( no more than roundoff-errors )
-You can use registers R56 to R64 to write the subroutine.

-Simply press R/S to continue with the same h & N or modify registers R05 & R06 if need be.

4°)  System of n Differential Equations

-We have to solve the system:

y"1 = f1(x,y1, ... ,yn)
y"2 = f2(x,y1, ... ,yn)
..................                                        with the initial values:     yi(x0)  &  y'i(x0)      i = 1 , ... , n

y"n = fn(x,y1, ... ,yn)

Data Registers:           •  R00 = subroutine name

( Registers R00 thru R03 , R20 thru R20+2n and R21+17n thru R136+17n are to be initialized before executing "NRKN10" )

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

•  R20 = x0
•  R21 = y1(x0)                                   •  R21+n = y'1(x0)                                   R21+2n thru R20+17n: temp
•  R22 = y2(x0)                                   •  R22+n = y'2(x0)
...............         ( initial values )           .......................

•  R20+n = yn(x0)                               •  R20+2n = y'n(x0)

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

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

•  R98+17n = a13,1  .................................    •  R109+17n = a13,12       •  R110+17n = b17

•  R111+17n = c1 ........................................................................       •  R123+17n = c17
•  R124+17n = c'1 .......................................................................       •  R136+17n = c'17

Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+2n, ... , R20+3n  respectively
with  x, y1, ... , yn in regiters R20, R21, ... , R20+n

-Line 162 is a three-byte GTO 01

 01 LBL "NRKN10"  02 RCL 03  03 STO 11  04 RCL 01  05 .1  06 %  07 +  08 STO 05  09 FRC  10 21.02  11 +  12 STO 04  13 RCL 05  14 +  15 STO 10  16 LASTX  17 ST+ X  18 +  19 STO 06  20 GTO 01  21 LBL 02  22 XEQ IND 00  23 RCL 04  24 RCL 05  25 ST+ X  26 +  27 RCL 07  28 LBL 00 29 RCL IND Y  30 RCL 02  31 *  32 STO IND Y  33 ISG Y  34 RDN  35 ISG Y  36 GTO 00  37 RTN  38 LBL 03  39 RCL IND Y  40 STO IND Y        41 ISG Y  42 RDN  43 ISG Y  44 GTO 03  45 RTN  46 LBL 07  47 RCL IND 07  48 RCL IND 08  49 *  50 RCL 02  51 *  52 ST+ IND 04  53 ISG 04  54 CLX  55 ISG 07  56 GTO 07 57 RCL 01  58 ST- 04  59 RCL 05  60 FRC  61 ST+ 07  62 ISG 08  63 GTO 07  64 RTN  65 LBL 01  66 RCL 05  67 RCL 06  68 +  69 STO 07              70 XEQ 02  71 RCL 04  72 RCL 06  73 XEQ 03  74 RCL 01  75 17  76 *  77 21  78 +  79 .1  80 %  81 +  82 STO 08  83 3.014  84 STO 09 85 RCL 20  86 STO 19  87 LBL 04  88 XEQ 07  89 RCL 02  90 RCL IND 08  91 *  92 ST+ 20  93 SIGN  94 LBL 06  95 RCL IND 10  96 LASTX  97 *  98 ST+ IND 04       99 ISG 04 100 CLX 101 ISG 10 102 GTO 06 103 RCL 01 104 ST- 04 105 ST- 10 106 XEQ 02 107 RCL 06 108 RCL 04 109 XEQ 03 110 RCL 19 111 STO 20 112 RCL 09 113 INT 114  E3 115 / 116 ISG X 117 ST+ 08 118 RCL 05 119 RCL 06 120 + 121 STO 07 122 ISG 09 123 GTO 04 124 RCL 02 125 ST+ 20 126 LBL 05 127 RCL IND 10    128 RCL 02 129 * 130 ST+ IND 04 131 ISG 10 132 CLX 133 ISG 04 134 GTO 05 135 RCL 01 136 ST- 04 137 ST- 10 138 XEQ 07 139 .013 140 ST+ 08 141 RCL 05 142 RCL 06 143 + 144 STO 07 145 LBL 08 146 RCL IND 07 147 RCL IND 08 148 * 149 ST+ IND 10 150 ISG 10 151 CLX 152 ISG 07 153 GTO 08            154 RCL 01 155 ST- 10 156 RCL 05 157 FRC 158 ST+ 07 159 ISG 08 160 GTO 08 161 DSE 11 162 GTO 01 163 RCL 23 164 RCL 22 165 RCL 21 166 RCL 20 167 END

( 273 bytes / SIZE 137+17.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:

d2y1/dx2 = - y1y2y3
d2y2/dx2 =  x ( y1+y2-y3 )
d2y3/dx2 =  xy1- y2y3

with the initial values:   y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2    &   y'1(0) = 1 ; y'2(0) = 1 ; y'3(0) = 1

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 R27 ; R28 ; R29

-For instance, the following subroutine:

 01  LBL "T"  02  RCL 21  03  RCL 22  04  RCL 23  05  *  06  *  07  CHS  08  STO 27  09  RCL 21  10  RCL 22  11  +  12  RCL 23  13  -  14  RCL 20  15  *  16  STO 28  17  RCL 20  18  RCL 21  19  *  20  RCL 22  21  RCL 23  22  *  23  -  24  STO29  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 116 constants into R21+17n to R136+17n
-Execute for instance the following routine:

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

-Then the initial values:

x  = 0   STO 20     and

y1 = 1   STO 21     y'1 =  1  STO 24
y2 = 1   STO 22     y'2 =  1  STO 25
y3 = 2   STO 23     y'3 =  1  STO 26

-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 "NRKN10"  >>>>        x   =  1                                                      ---Execution time = 4s---            with V41 in turbo mode

RDN     y1(1) =  0.439524097 = R21                     y'1(1) = -2.101122877 = R24
RDN     y2(1) =  2.070940658 = R22       and         y'2(1) =  1.269596947 = R25
RDN     y3(1) =  1.744524962 = R23                      y'3(1) = -1.704234764 = R26

Notes:

-With a real HP41, execution time is about 40 minutes, but with free42, the results are returned almost instantaneously !
-Since SIZE = 137+17n, "NRKN10" can solve at most a system of n = 10 differential equations with SIZE 307.
-But with free42, n may reach 50 with SIZE 990 ( beyond register 999, the control numbers would not work )

Reference:

[1]  Philip Sharp - http://www.math.auckland.ac.nz/~sharp/rkn1.dat

-The programs listed in this page employ formulae without error-estimate, so you'll have to calculate again
the solutions with a smaller stepsize - say  h/2 - to get an idea of the precision.
-But in reference [1], you will find the coefficients of pairs of embedded Runge-Kutta-Nystrom formulas
of various orders to control the stepsize more directly.