The Museum of HP Calculators

# Numerov's Method 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°) Numerov's method

a)  1 differential equation      y" = f(x,y)
b)  2 differential equations    y" = f(x,y,z)  ;  z" = g(x,y,z)
c)  n differential equations    y"1 = f1(x;y1,.....,yn)  , ....... , y"n = fn(x;y1,.....,yn)     ( X-Functions Module required )

2°) A formula of order 7

a)  1 differential equation      y" = f(x,y)
b)  2 differential equations    y" = f(x,y,z)  ;  z" = g(x,y,z)
c)  n differential equations    y"1 = f1(x;y1,.....,yn)  , ....... , y"n = fn(x;y1,.....,yn)      ( X-Functions Module required )

-Numerov's method has been devised to solve second order differential equations where the first derivative y' does not appear explicitly.
-The formula is:  yn+1 = 2.yn - yn-1 + (h2/12).( fn+1 + 10.fn + fn-1 )    where h = xn - xn-1 is the stepsize and  fk = f(xk,yk)  ;  yn-j = y(xn-j.h)

-It requires 2 starting values  y-1 and y0 instead of  y0 and y'0 .
-The method duplicates the Taylor series up to the term in h5

-Since it is an implicit formula ( yn+1 appears on both sides ) , we use an iterative method at each step.

1°) Numerov's Method

a) 1 differential equation    y" = f(x,y)

Data Registers:           •  R00 =  f name                                       ( Registers R00 thru R05 are to be initialized before executing "NUM1" )

•  R01 = x0
•  R02 = y0   •  R03 = y-1 = y(x0 - h )    •  R04 = h = stepsize    •  R05 = N = number of steps      R06 to R10: temp
Flags: /
Subroutine:   A program which computes  f(x,y)  assuming  x is in X-register and y is in Y-register upon entry

01  LBL "NUM1"
02  RCL 05
03  STO 06
04  RCL 03
05  RCL 01
06  RCL 04
07  -
08  XEQ IND 00
09  STO 07
10  RCL 02
11  STO 10
12  RCL 01
13  XEQ IND 00
14  STO 08
15  LBL 01
16  RCL 10
17  RCL 01
18  RCL 04
19  +
20  XEQ IND 00
21  STO 09
22  RCL 08
23  10
24  *
25  +
26  RCL 07
27  +
28  RCL 04
29  X^2
30  *
31  12
32  /
33  RCL 02
34  ST+ X
35  RCL 03
36  -
37  +
38  ENTER^
39  X<> 10
40  VIEW X             this line is only useful to display the successive approximations, but it's not necessary for the calculations.
41  X#Y?
42  GTO 01
43  X<> 02
44  STO 03
45  RCL 04
46  ST+ 01
47  RCL 09
48  X<> 08
49  STO 07
50  DSE 06
51  GTO 01
52  RCL 02
53  RCL 01
54  CLD
55  END

( 78 bytes / SIZE 011 )

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

Example:     y" = ( x2 - 1 ).y   ;   y(0) = 1 &  y(-0.1) = 0.995012479

-Let's evaluate y(1)

-First, we load the required subroutine in main memory, for instance:

LBL "T"     any global label, maximum 6 characters
X^2
1
-
*
RTN

"T"  ASTO 00

0   STO 01
1   STO 02      0.995012479   STO 03
0.1  STO 04                       10   STO 05     ( here,  h = 0.1 and  N = 10 )

XEQ "NUM1"  gives       x   = 1                           ( in 53 seconds ( or 48 seconds if you delete line 40 ) )
X<>Y   y(1) = 0.606528753

-To continue with the same N-value, simply press R/S , it yields  y(2) = 0.135332761

-The solution was actually  y(x) = exp ( -x2/2 )   so  y(1) = 0.606530660  and  y(2) = 0.135335283  the error is of the order of  3 E-6

Notes:

-In a trajectory problem , we can use tabulated values to begin the calculations.
-If we don't have 2 starting values, the second may be found by a Runge-Kutta method, provided you know y'0 .

-The loop stops when 2 consecutive approximations are equal ( line 41 ).
-This might lead to an infinite loop but practically, the risk is tiny because of the factor h2 in the formula - provided h remains "reasonably small" however.

b) 2 differential equations    y" = f(x,y,z)  ;   z" = g(x,y,z)

-Of course, the method may be generalized to systems of differential equations:

Data Registers:           •  R00 =  subroutine name                       ( Registers R00 thru R07 are to be initialized before executing "NUM2" )

•  R01 = x0
•  R02 = y0   •  R04 = h = stepsize                 •  R06 = y-1 = y(x0 - h )
•  R03 = z0   •  R05 = N = number of steps    •  R07 = z-1 = z(x0 - h )              R08 to R16: temp
Flags: /
Subroutine:   A program which computes  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 "NUM2"
02  RCL 05
03  STO 16
04  RCL 07
05  RCL 06
06  RCL 01
07  RCL 04
08  -
09  XEQ IND 00
10  STO 13
11  X<>Y
12  STO 10
13  RCL 03
14  STO 09
15  RCL 02
16  STO 08
17  RCL 01
18  XEQ IND 00
19  STO 14
20  X<>Y
21  STO 11
22  LBL 01
23  RCL 09
24  RCL 08
25  RCL 01
26  RCL 04
27  +
28  XEQ IND 00
29  STO 15
30  RCL 14
31  10
32  *
33  +
34  RCL 13
35  +
36  X<>Y
37  STO 12
38  RCL 11
39  10
40  *
41  +
42  RCL 10
43  +
44  RCL 04
45  X^2
46  12
47  /
48  ST* Z
49  *
50  RCL 06
51  -
52  RCL 02
53  ST+ X
54  +
55  ENTER^
56  X<> 08
57  -
58  ABS
59  X<>Y
60  RCL07
61  -
62  RCL 03
63  ST+ X
64  +
65  ENTER^
66  X<> 09
67  -
68  ABS
69  +
70  X#0?
71  GTO 01
72  RCL 08
73  X<> 02
74  STO 06
75  RCL 09
76  X<> 03
77  STO 07
78  RCL 12
79  X<> 11
80  STO 10
81  RCL 15
82  X<> 14
83  STO 13
84  RCL 04
85  ST+ 01
86  DSE 16
87  GTO 01
88  RCL 03
89  RCL 02
90  RCL 01
91  END

( 120 bytes / SIZE 017 )

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

Example:      y" = ( x - 2 ).z       y(1) = 0.367879441      y(0.9) = 0.365912694
z" = y/x                 z(1) = 0.367879441      z(0.9) = 0.406569660

-Evaluate  y(2) & z(2)

LBL "T"
ST/ Y
2
-
ST* Z
RDN
RTN

"T"  ASTO 00
1    STO 01      0.367879441   STO 02    STO 03                  h = 0.1  STO 04      ,       N = 10  STO 05
0.365912694   STO 06
0.406569660   STO 07

XEQ "NUM2"  >>>>       x = 2                        ( in 81 seconds )
RDN  y(2)  = 0.270670254
RDN  z(2)  = 0.135335322

-In fact, y = x exp(-x)  and  z = exp(-x) , so the exact results are  y(2)  = 0.270670566  &  z(2)  = 0.135335283

-Press R/S to continue the calculations ( after changing N in R05 if needed )

c)  n differential equations    y"1 = f1(x;y1,.....,yn)  , ....... , y"n = fn(x;y1,.....,yn)

Data Registers:           •  R00 =  subroutine name       ( Registers R00 thru R02 and R09 thru R2n+10 are to be initialized before executing "NUM" )

•  R01 = N = number of steps    •  R02 = h = stepsize            R03 to R08: temp

•  R09 = n = number  of equations
•  R10 = x0
•  R11 = y1(x0)               •  R11+n = y1(x0-h)
•  R12 = y2(x0)               •  R12+n = y2(x0-h)                         R11 thru R10+2n  contain
.....................                 ..........................                            the 2 starting values

•  R10+n = yn(x0)           •  R10+2n = yn(x0-h)

-During the calculations,    R11+n to R10+2n = fi(1)          the n values of the functions for  x =  x0+h
R11+2n to R10+3n = fi(0)        -------------------------------  x =  x0
R11+3n to R10+4n = fi(-1)      -------------------------------  x =  x0-h
R11+4n to R10+5n = yi(0)       --------------- yi(x0)
R11+5n to R10+6n = yi(-1)      --------------- yi(x0-h)

Flags: /
Subroutine:   A program which computes and stores   f1(x;y1,.....,yn)  , ....... ,  fn(x;y1,.....,yn) into  R11+n , .......... , R10+2n   respectively
with  x ; y1,.....,yn  in R10 ; R11 , ......... , R10+n   respectively

01  LBL "NUM"
02  RCL 09
03  4
04  *
05  11.011
06  STO 07
07  INT
08  +
09  STO 03
10   E3
11  /
12  RCL 07
13  INT
14  +
15  RCL 09
16   E6
17  /
18  STO 04
19  ST+ X
20  +
21  REGMOVE
22  XEQ IND 00
23  RCL 09
24  .1
25  %
26  ST+ X
27  +
28  RCL 07
29  +
30  RCL 04
31  +
32  STO 05
33  REGMOVE
34  RCL 03
35  RCL 07
36  FRC
37  +
38  RCL 04
39  +
40  STO 06
41  RCL 09
42  +
43  REGMOVE
44  RCL 02
45  ST- 10
46  XEQ IND 00
47  RCL 05
48  RCL 09
49   E3
50  /
51  +
52  REGMOVE
53  RCL 06
54  REGMOVE
55  RCL 02
56  ST+ 10                        Lines 01 to 56 are only executed the first time in order to initialize some registers.
57  LBL 01
58  RCL 01
59  STO 08
60  LBL 02
61  RCL 02
62  ST+ 10
63  LBL 03
64  XEQ IND 00
65  RCL 09
66  RCL 09
67  RCL 09
68  10.01
69  +
70  STO 07
71  INT
72  +
73  STO 06
74  +
75  STO 05
76  +
77  STO 04
78  +
79  STO 03
80  CLX
81  LBL 04
82  RCL IND 04
83  RCL IND 05
84  10
85  *
86  +
87  RCL IND 06
88  +
89  RCL 02
90  X^2
91  *
92  12
93  /
94  RCL IND 03
95  ST+ X
96  +
97  RCL 03
98  RCL 09
99  +
100  RDN
101  RCL IND T
102  -
103  ENTER^
104  X<> IND 07
105  -
106  ABS
107  +
108  DSE 03
109  DSE 04
110  DSE 05
111  DSE 06
112  DSE 07
113  GTO 04
114  X#0?
115  GTO 03                                the loop stops when 2 consecutive approximations are equal
116  RCL 05
117  1
118  +
119  .1
120  %
121  +
122  RCL 09
123  ST- Y
124   E6
125  /
126  STO 07
127  4
128  *
129  +
130  REGMOVE
131  RCL 03
132  1
133  +
134   E3
135  /
136  RCL 07
137  +
138  11
139  +
140  REGMOVE
141  DSE 08
142  GTO 02
143  RCL 13
144  RCL 12
145  RCL 11
146  RCL 10
147  RTN
148  GTO 01
149  END

( 210 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:     We want to study the motion of a planet around a point sun ( the mass of the planet is neglected ).
Here, the variable is t = time ( in days ), and  y1 , y2 , y3  are denoted  x , y , z  ( the coordinates are expressed in Astronomical Units )

x" = -k2 x/(x2+y2+z2)3/2
-The equations are      y" = -k2 y/(x2+y2+z2)3/2          where  k = 0.01720209895  is Gauss' constant
z" = -k2 z/(x2+y2+z2)3/2

-At  t = 0   x =  0.092     at  t = -1   x =  0.070      Find the position of the planet at  t = 2 days
y = -0.445                     y = -0.451
z = -0.045                     z = -0.043

-We store  -k2  in register R30        17.20209895 E-3   X^2  CHS  STO 30
-we load the following subroutine in main memory:

LBL "T"      X^2          RCL 13        SQRT        STO 14      RCL 29       RCL 30      STO 16
RCL 30      RCL 12    X^2              *                RCL 12       /                  *                RTN
RCL 11      X^2          +                  STO 29      RCL 30      STO 15       RCL 29
ST* Y        +               ENTER^      /                 *                 RCL 13       /

"T"  ASTO 00     number of steps = 2  STO 01     stepsize h = 1  STO 02    there are 3 equations:   3  STO 09    then the initial values:

0       STO 10
0.092   STO 11      0.070    STO 14
-0.445   STO 12     -0.451    STO 15
-0.045   STO 13     -0.043    STO 16

XEQ "NUM"  >>>>   t =   2                                      ( execution time = 56 seconds )
RDN   x =   0.135070  =  R11
RDN   y = -0.428856  =  R12
RDN   z = -0.048573  =  R13

N.B:    To continue, press R/S or XEQ 01 ( after changing N in register R01 if needed )
but not   XEQ "NUM"  because the second new initial values are not in registers  R14 thru R16 but in R26-R27-R28

-With the same N-value ( N = 2 ), you should find   t = 4 ,  x = 0.176408 , y = -0.407227 , z = -0.051524

-For a planet like Mercury with h =  1   day, the error is of the order of   2.7 10 -5  AU  after  88 days
---- h = 0.5 day, --------------------------   1.6 10 -6  AU  -------------

-When h is divided by 2, the errors are approximately divided by 16 = 24

2°) A formula of order 7

-We can use more than 2 starting values in order to reach a better accuracy.
-The following formula duplicates the Taylor's series up to the term in h7

yn+1 = yn + yn-2 - yn-3  + (h2/240).( 17. fn+1 + 232.fn + 222.fn-1 + 232.fn-2 + 17.fn-3 )       [  h = the stepsize ,  fk = f(xk,yk)  ,  yn-j = y(xn-j.h)  ]

-4 values are needed to initialize the algorithm:  y0 , y-1 , y-2 , y-3

a) 1 differential equation    y" = f(x,y)

Data Registers:           •  R00 =  f name                                       ( Registers R00 thru R07 are to be initialized before executing "7NUM1" )

•  R01 = x0
•  R02 = y0   •  R03 = y-1 = y(x0 - h )       •  R04 = h = stepsize    •  R05 = N = number of steps
•  R06 = y-2 = y(x0 - 2.h )
•  R07 = y-3 = y(x0 - 3.h )                                            R08 to R14: temp
Flags: /
Subroutine:   A program which computes  f(x,y)  assuming  x is in X-register and y is in Y-register upon entry

01  LBL "7NUM1"
02  RCL 05
03  STO 08
04  RCL 07
05  RCL 01
06  RCL 04
07  3
08  *
09  -
10  XEQ IND 00
11  STO 09
12  RCL 06
13  RCL 01
14  RCL 04
15  ST+ X
16  -
17  XEQ IND 00
18  STO 10
19  RCL 03
20  RCL 01
21  RCL 04
22  -
23  XEQ IND 00
24  STO 11
25  RCL 02
26  STO 14
27  RCL 01
28  XEQ IND 00
29  STO 12
30  LBL 01
31  RCL 14
32  RCL 01
33  RCL 04
34  +
35  XEQ IND 00
36  STO 13
37  RCL 09
38  +
39  17
40  *
41  RCL 10
42  RCL 12
43  +
44  232
45  *
46  +
47  RCL 11
48  222
49  *
50  +
51  240
52  /
53  RCL 04
54  X^2
55  *
56  RCL 07
57  -
58  RCL 06
59  +
60  RCL 02
61  +
62  ENTER^
63  X<> 14
64  X#Y?
65  GTO 01
66  X<> 02
67  X<> 03
68  X<> 06
69  STO 07
70  RCL 13
71  X<> 12
72  X<> 11
73  X<> 10
74  STO 09
75  RCL 04
76  ST+ 01
77  DSE 08
78  GTO 01
79  RCL 02
80  RCL 01
81  END

( 115 bytes / SIZE 015 )

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

Example:     y" = ( x2 - 1 ).y   Knowing   y(0) = 1 ;  y(-0.1) = 0.995012479 ;  y(-0.2) = 0.980198673 ;  y(-0.3) = 0.955997482

-Evaluate y(1)

LBL "T"
X^2
1
-
*
RTN

"T"  ASTO 00

0   STO 01
1   STO 02       0.995012479   STO 03           h =  0.1  STO 04  ,  N = 10   STO 05
0.980198673   STO 06
0.955997482   STO 07

XEQ "7NUM1"  gives       x   = 1                           ( in 70 seconds )
X<>Y   y(1) = 0.606530689

-To continue with the same N-value, simply press R/S , it yields  y(2) = 0.135335319

-The errors are of the order of  3 E-8

b) 2 differential equations    y" = f(x,y,z)  ;   z" = g(x,y,z)

Data Registers:           •  R00 =  subroutine name                       ( Registers R00 thru R11 are to be initialized before executing "7NUM2" )

•  R01 = x0
•  R02 = y0   •  R04 = h = stepsize                 •  R06 = y-1 = y(x0 - h )          •  R09 = z-1 = z(x0 - h )
•  R03 = z0   •  R05 = N = number of steps    •  R07 = y-2 = y(x0 - 2h )        •  R10 = z-2 = z(x0 - 2h )           R12 to R24: temp
•  R08 = y-3 = y(x0 - 3h )        •  R11 = z-3 = z(x0 - 3h )
Flags: /
Subroutine:   A program which computes  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 "7NUM2"
02  RCL 05
03  STO 24
04  RCL 11
05  RCL 08
06  RCL 04
07  3
08  *
09  RCL 01
10  X<>Y
11  -
12  XEQ IND 00
13  STO 17
14  X<>Y
15  STO 12
16  RCL 10
17  RCL 07
18  RCL 01
19  RCL 04
20  ST+ X
21  -
22  XEQ IND 00
23  STO 18
24  X<>Y
25  STO 13
26  RCL 09
27  RCL 06
28  RCL 01
29  RCL 04
30  -
31  XEQ IND 00
32  STO 19
33  X<>Y
34  STO 14
35  RCL 03
36  STO 23
37  RCL 02
38  STO 22
39  RCL 01
40  XEQ IND 00
41  STO 20
42  X<>Y
43  STO 15
44  LBL 01
45  RCL 23
46  RCL 22
47  RCL 01
48  RCL 04
49  +
50  XEQ IND 00
51  STO 21
52  RCL 17
53  +
54  17
55  *
56  RCL 18
57  RCL 20
58  +
59  232
60  *
61  +
62  RCL 19
63  222
64  *
65  +
66  X<>Y
67  STO 16
68  RCL 12
69  +
70  17
71  *
72  RCL 13
73  RCL 15
74  +
75  232
76  *
77  +
78  RCL 14
79  222
80  *
81  +
82  RCL 04
83  X^2
84  240
85  /
86  ST* Z
87  *
88  RCL 08
89  -
90  RCL 07
91  +
92  RCL 02
93  +
94  ENTER^
95  X<> 22
96  -
97  ABS
98  X<>Y
99  RCL 11
100  -
101  RCL 10
102  +
103  RCL 03
104  +
105  ENTER^
106  X<> 23
107  -
108  ABS
109  +
110  X#0?
111  GTO 01
112  RCL 22
113  X<> 02
114  X<> 06
115  X<> 07
116  STO 08
117  RCL 23
118  X<> 03
119  X<> 09
120  X<> 10
121  STO 11
122  RCL 16
123  X<> 15
124  X<> 14
125  X<> 13
126  STO 12
127  RCL 21
128  X<> 20
129  X<> 19
130  X<> 18
131  STO 17
132  RCL 04
133  ST+ 01
134  DSE 24
135  GTO 01           a three-byte GTO
136  RCL 03
137  RCL 02
138  RCL 01
139  END

( 207 bytes / SIZE 025 )

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

Example:      y" = ( x - 2 ).z        y(1) = 0.367879441      y(0.9) = 0.365912694     y(0.8) = 0.359463171    y(0.7) = 0.347609713
z" = y/x                  z(1) = 0.367879441      z(0.9) = 0.406569660     z(0.8) = 0.449328964    z(0.7) =  0.496585304

-Evaluate  y(2) & z(2)

LBL "T"
ST/ Y
2
-
ST* Z
RDN
RTN

"T"  ASTO 00
1    STO 01      0.367879441  STO 02    STO 03                  h = 0.1  STO 04       ,      N = 10  STO 05

0.365912694   STO 06           0.406569660   STO 09
0.359463171   STO 07           0.449328964   STO 10
0.347609713   STO 08           0.496585304   STO 11

XEQ "7NUM2"  >>>>     x = 2                        ( in 125 seconds )
RDN  y(2)  = 0.270670563
RDN  z(2)  = 0.135335281

-The errors are of the order of  3 E-9

-Press R/S to continue the calculations ( after changing N in R05 if needed )

c)  n differential equations    y"1 = f1(x;y1,.....,yn)  , ....... , y"n = fn(x;y1,.....,yn)

Data Registers:           •  R00 =  subroutine name       ( Registers R00 thru R02 and R13 thru R4n+14 are to be initialized before executing "7NUM" )

•  R01 = N = number of steps    •  R02 = h = stepsize            R03 to R12: temp

•  R13 = n = number  of equations
•  R14 = x0
•  R15 = y1(x0)               •  R15+n = y1(x0-h)               •  R15+2n = y1(x0-2h)               •  R15+3n = y1(x0-3h)
•  R16 = y2(x0)               •  R16+n = y2(x0-h)               •  R16+2n = y2(x0-2h)               •  R16+3n = y2(x0-3h)      R15 thru R14+4n
.....................                 ..........................                   ..............................                                                                  contain

•  R14+n = yn(x0)           •  R14+2n = yn(x0-h)             •  R14+3n = yn(x0-2h)              •  R14+4n = yn(x0-3h)      the 4 starting values

-During the calculations,    R15+n to R14+2n = fi(1)          the n values of the functions for  x =  x0+h
R15+2n to R14+3n = fi(0)        -------------------------------  x =  x0
R15+3n to R14+4n = fi(-1)       -------------------------------  x =  x0-h
R15+4n to R14+5n = fi(-2)       -------------------------------  x =  x0-2h
R15+5n to R14+6n = fi(-3)       -------------------------------  x =  x0-3h
R15+6n to R14+7n = yi(0)       --------------- yi(x0)
R15+7n to R14+8n = yi(-1)      --------------- yi(x0-h)
R15+8n to R14+9n = yi(-2)      --------------- yi(x0-2h)
R15+9n to R14+10n = yi(-3)    --------------- yi(x0-3h)

Flags: /
Subroutine:   A program which computes and stores   f1(x;y1,.....,yn)  , ....... ,  fn(x;y1,.....,yn) into  R15+n , .......... , R14+2n   respectively
with  x ; y1,.....,yn  in R14 ; R15 , ......... , R14+n   respectively

01  LBL "7NUM"
02  RCL 14
03  STO 10
04  RCL 13
05   E6
06  /
07  STO 03
08  4
09  *
10  RCL 13
11   E3
12  /
13  STO 04
14  6
15  *
16  +
17  15.015
18  STO 05
19  +
20  REGMOVE
21  LASTX
22  RCL 03
23  +
24  RCL 04
25  +
26  RCL 13
27  +
28  STO 06
29  RCL 03
30  RCL 05
31  +
32  RCL 13
33  6
34  *
35  +
36  STO 07
37  STO 08
38  4
39  STO 09
40  LBL 10
41  XEQ IND 00
42  RCL 04
43  RCL 06
44  +
45  STO 06
46  REGMOVE
47  DSE 09
48  X=0?
49  GTO 00
50  RCL 07
51  RCL 13
52  +
53  STO 07
54  REGMOVE
55  RCL 02
56  ST- 14
57  GTO 10
58  LBL 00
59  RCL 08
60  REGMOVE
61  RCL 10
62  STO 14                        Lines 01 to 62 are only executed the first time in order to initialize some registers.
63  LBL 01
64  RCL 01
65  STO 12
66  LBL 02
67  RCL 02
68  ST+ 14
69  LBL 03
70  XEQ IND 00
71  RCL 13
72  RCL 13
73  RCL 13
74  14.014
75  +
76  STO 03
77  INT
78  +
79  STO 04
80  +
81  STO 05
82  +
83  STO 06
84  +
85  STO 07
86  +
87  STO 08
88  +
89  STO 09
90  +
91  +
92  STO 10
93  +
94  STO 11
95  CLX
96  LBL 04
97  RCL IND 04
98  RCL IND 08
99  +
100  17
101  *
102  RCL IND 05
103  RCL IND 07
104  +
105  232
106  *
107  +
108  RCL IND 06
109  222
110  *
111  +
112  RCL 02
113  X^2
114  *
115  240
116  /
117  RCL IND 09
118  +
119  RCL IND 10
120  +
121  RCL IND 11
122  -
123  ENTER^
124  X<> IND 03
125  -
126  ABS
127  +
128  DSE 11
129  DSE 10
130  DSE 09
131  DSE 08
132  DSE 07
133  DSE 06
134  DSE 05
135  DSE 04
136  DSE 03
137  GTO 04
138  X#0?
139  GTO 03
140  RCL 05
141  1
142  +
143  .1
144  %
145  +
146  RCL 13
147  ST- Y
148   E6
149  /
150  STO 07
151  8
152  *
153  +
154  REGMOVE
155  RCL 09
156  1
157  +
158   E3
159  /
160  RCL 07
161  +
162  15
163  +
164  REGMOVE
165  DSE 12
166  GTO 02                   a three-byte GTO
167  RCL 17
168  RCL 16
169  RCL 15
170  RCL 14
171  RTN
172  GTO 01                   a three-byte GTO
173  END

( 246 bytes / SIZE 10n+15 )

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

Example:     We study again the motion of a planet around a point sun ( the mass of the planet is neglected ).
The variable is t = time ( in days ), and  y1 , y2 , y3  are denoted  x , y , z  ( the coordinates are expressed in Astronomical Units )

x" = -k2 x/(x2+y2+z2)3/2
-The equations are      y" = -k2 y/(x2+y2+z2)3/2          where  k = 0.01720209895  is Gauss' constant
z" = -k2 z/(x2+y2+z2)3/2

-At  t = 0     x =  0.293510249       at  t = -1     x =  0.301200207      at t = -2     x = 0.305864609       at  t = -3     x = 0.307427938
y =  0.091967806                          y =  0.061830391                       y =  0.031072548                         y = 0
z =  0.040946705                          z =  0.027528664                       z =  0.013834390                          z = 0

Find the position of the planet at  t = 4 days

-We store  -k2  in register R45        17.20209895 E-3   X^2  CHS  STO 45
and we load the following subroutine in main memory:

LBL "T"      X^2          RCL 17        SQRT        STO 18      RCL 46       RCL 45      STO 20
RCL 45      RCL 16    X^2              *                RCL 16       /                  *                RTN
RCL 15      X^2          +                  STO 46      RCL 45      STO 19       RCL 46
ST* Y        +               ENTER^      /                 *                 RCL 17       /

"T"  ASTO 00     number of steps = 4  STO 01     stepsize h = 1  STO 02    there are 3 equations:   3  STO 13    then the initial values:

0            STO 14
0.293510249   STO 15      0.301200207    STO 18    0.305864609     STO 21   0.307427938   STO 24
0.091967806   STO 16      0.061830391    STO 19    0.031072548     STO 22             0            STO 25
0.040946705   STO 17      0.027528664    STO 20    0.013834390     STO 23             0            STO 26

XEQ "7NUM"  >>>   t =  4                     =  R14                    ( execution time = 2m32s )
RDN   x =  0.235500989  =  R15
RDN   y =  0.200940664  =  R16
RDN   z =  0.089464547  =  R17

N.B:    To continue, press R/S or XEQ 01 ( after changing N in register R01 if needed )  but not   XEQ "7NUM"

-For a planet like Mercury, with h =  1   day, the error is of the order of   3.6 10 -7  AU  after  88 days
---- h = 0.5 day, --------------------------   5.8 10 -9  AU  -------------

-When h is divided by 2, the errors are approximately divided by 64 = 26

-In fact, the error - after 88 days with  h = 0.5 day - is 5.8 10 -9  AU  on an HP-48,
but it reaches 1.6 10 -8 AU on an HP-41 because of roundoff errors: the HP-41 works with 10 digits only!

Reference:

[1]  Francis Scheid   "Numerical Analysis"       McGraw-Hill   ISBN 07-055197-9