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