hp41programs

Numerov

Numerov's Method for the HP-41


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

-Line 40 is only useful to display the successive approximations, but it's not necessary for the calculations.
 
 

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 
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)        the n values of the functions for  x =  x0
                R11+3n to R10+4n = fi(-1)      the n values of the functions for  x =  x0-h
                R11+4n to R10+5n = yi(0)       the n values of  yi(x0)
                R11+5n to R10+6n = yi(-1)      the n values of  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

-Lines 01 to 56 are only executed the first time in order to initialize some registers.
 
 

  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    
  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  
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
      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
       y = -0.445
       z = -0.045

-At  t = -1

       x =  0.070                               Find the position of the planet at  t = 2 days
       y = -0.451
       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
    with h = 0.5 day, the error is of the order of    1.6 10 -6  AU  after  88 days

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

-Line 135 is a three-byte  GTO 01
 
 

  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 
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 =
                                        .....................                 ..........................                   ..............................

                                      •  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)        the n values of the functions for  x =  x0
            R15+3n to R14+4n = fi(-1)       the n values of the functions for  x =  x0-h
            R15+4n to R14+5n = fi(-2)       the n values of the functions for  x =  x0-2h
            R15+5n to R14+6n = fi(-3)       the n values of the functions for  x =  x0-3h
            R15+6n to R14+7n = yi(0)       the n values of  yi(x0)
            R15+7n to R14+8n = yi(-1)      the n values of  yi(x0-h)
            R15+8n to R14+9n = yi(-2)      the n values of  yi(x0-2h)
            R15+9n to R14+10n = yi(-3)    the n values of  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

-Lines 01 to 62 are only executed the first time in order to initialize some registers.
-Lines 166 & 172 are three-byte  GTOs
 
 

  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  
  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 
167  RCL 17
168  RCL 16
169  RCL 15
170  RCL 14
171  RTN
172  GTO 01 
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
      y" = -k2 y/(x2+y2+z2)3/2          where  k = 0.01720209895  is Gauss' constant
      z" = -k2 z/(x2+y2+z2)3/2

                    t = 0                            t = -1                            t = -2                            t = -3

       x =  0.293510249          x =  0.301200207          x = 0.305864609           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
   with h = 0.5 day, the error is of the order of   5.8 10 -9  AU  after  88 days

-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