# hp41programs

Diffeq Systems of Differential Equations for the HP-41

Overview

1°) Systems of first-order Differential Equations

a) The "Classical" 4th-order Runge-Kutta Method
b) A Runge-Kutta Method of order 6
c) A Runge-Kutta Method of order 8
d) Runge-Kutta-Fehlberg 4-5 Method
e) Bulirsch-Stoer Method

2°) Systems of second-order Conservative Equations

a)  Systems of 2 Equations
b)  Systems of 3 Equations
c)  Systems of n Equations

-The programs listed in paragraph 1°)  solve a system of n first-order differential equations:

dy1/dx = f1(x,y1, ... ,yn)
dy2/dx = f2(x,y1, ... ,yn)
..................                                        with the initial values:     yi(x0)         i = 1 , ... , n

dyn/dx = fn(x,y1, ... ,yn)

-In paragraph 2°) the HP-41 solves systems of second-order equations where the first derivatives  y'i(x) do not appear explicitly:

d2y1/dx2 = f1(x,y1, ... ,yn)
d2y2/dx2 = f2(x,y1, ... ,yn)
..................                                        with the initial values:     yi(x0)   ,    y'i(x0)       i = 1 , ... , n

d2yn/dx2 = fn(x,y1, ... ,yn)

-All these programs use 2 instructions from the X-functions module:  REGMOVE  and  REGSWAP - Except in §2°)a)b)

1°) Systems of first-order Differential Equations

a) The "Classical" 4th-order Runge-Kutta Method

Data Registers:           •  R00 = subroutine name       ( Registers R00 thru R03 and R10 thru R10+n are to be initialized before executing "RK4" )

•  R01 = n  = number of equations = number of functions         R04 thru R09: temp
•  R02 = h  = stepsize
•  R03 = N = number of steps

•  R10 = x0
•  R11 = y1(x0)                                                                           R11+n thru R10+4n: temp
•  R12 =  y2(x0)
...............         ( initial values )

•  R10+n = yn(x0)
Flags:  F06-F07
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R11+n, ... , R10+2n  respectively
with  x, y1, ... , yn in regiters R10, R11, ... , R10+n

-If you don't have an HP-41 CX, line 24 can be replaced by the 5 lines:

0
LBL 00
STO IND Y
ISG Y
GTO 00

 01  LBL "RK4"  02  RCL 03 03  STO 09 04  RCL 01 05  RCL 01 06  RCL 01 07  10.01 08  + 09  STO 04 10  INT 11  + 12  STO 05 13  + 14  STO 06 15  + 16  STO 07 17  RCL 05 18  RCL 06 19  1 20  + 21   E3 22  / 23  + 24  CLRGX 25  11 26  LASTX 27  + 28  RCL 01 29   E6 30  / 31  + 32  STO 08 33  GTO 03 34  LBL 01 35  XEQ IND 00 36  LBL 02 37  RCL IND 05 38  RCL 02 39  * 40  ST+ IND 06 41  FS? 06 42  ST+ IND 06 43  FC? 07 44  ST+ X 45  RCL IND 07 46  + 47  STO IND 04 48  DSE 07 49  DSE 06 50  DSE 05 51  DSE 04 52  GTO 02 53  RCL 01 54  ST+ 04 55  ST+ 05 56  ST+ 06 57  ST+ 07 58  RTN 59  LBL 03 60  RCL 08 61  REGMOVE 62  CF 06 63  SF 07 64  XEQ 01 65  SF 06 66  RCL 02 67  ST+ 10 68  XEQ 01 69  CF 07 70  XEQ 01 71  CF 06 72  RCL 02 73  ST+ 10 74  XEQ 01 75  RCL 08 76  REGSWAP 77  RCL 04 78  RCL 06 79  3 80  SIGN 81  LBL 04 82  CLX 83  X<> IND Y 84  LASTX 85  / 86  ST+ IND Z 87  DSE Y 88  DSE Z 89  GTO 04 90  DSE 09 91  GTO 03 92  RCL 13 93  RCL 12 94  RCL 11 95  RCL 10 96  END

( 155 bytes / SIZE 4n+11 )

 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 the system:

dy1/dx = - y1y2y3
dy2/dx =  x ( y1+y2-y3 )
dy3/dx =  xy1- y2y3

with the initial values:   y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2

x ; y1 ; y2 ; y3  will be in R10 ; R11 ; R12 ; R13 respectively, and we have to write a program to compute the 3 functions
and store the results in R14 ; R15 ; R16   ( just after the "initial-value registers" ).

-For instance, the following subroutine will do the job - let's call it "T"

01  LBL "T"
02  RCL 11
03  RCL 12
04  RCL 13
05  *
06  *
07  CHS
08  STO 14
09  RCL 11
10  RCL 12
11  +
12  RCL 13
13  -
14  RCL 10
15  *
16  STO 15
17  RCL 10
18  RCL 11
19  *
20  RCL 12
21  RCL 13
22  *
23  -
24  STO16
25  RTN

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

then the initial values:

x     in R10                      0 STO 10
y1   in R11                      1 STO 11
y2   in R12                      1 STO 12
y3   in R13                      2 STO 13

-Suppose we want to know   y1(1) , y2(1) ,  y3(1)  with a step size h = 0.1 and therefore N = 10

h/2 is stored in R02                0.05  STO 02
N is stored in R03                 10    STO 03   and  XEQ "RK4"

about 2mn24s later, the HP-41 returns:

x  =  1                        in X and R10
y1(1) =  0.2582093855  in Y and R11
y2(1) =  1.157619553    in Z and R12
y3(1) =  0.8421786516  in T and R13

-To continue with the same h and N, simply press R/S and you'll have y1(2) y2(2) y3(2)

-An estimation of the accuracy may be obtained by doing the calculation again with a smaller step size like h=0.05

It yields:       y1(1) = 0.2582079989    y2(1) = 1.157623732    y3(1) = 0.8421783424

Theoretically, the results tend to the exact solution as h tends to zero,
but naturally, roundoff errors ( and execution time ) will limit the accuracy attainable!

b) A Runge-Kutta Method of Order 6

Data Registers:           •  R00 = subroutine name       ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RK6" )

•  R01 = n  = number of equations = number of functions        R04 thru R13: temp        R14 thru R19: unused
•  R02 = h  = stepsize
•  R03 = N = number of steps

•  R20 = x0
•  R21 = y1(x0)                                                                           R21+n thru R20+8n: temp
•  R22 =  y2(x0)
...............         ( initial values )

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

-Line 301 is a three-byte  GTO 01

 01  LBL "RK6"   02  RCL 03   03  STO 04   04  21.021   05  RCL 01   06   E6   07  /   08  +   09  RCL 01   10  7   11  *   12  +   13  STO 05   14  RCL 01   15  RCL 01   16  RCL 01   17  20.02   18  +   19  STO 06   20  +   21  STO 07   22  +   23  STO 08   24  +   25  STO 09   26  +   27  STO 10   28  +   29  STO 11   30  +   31  STO 12   32  LBL 01   33  RCL 20   34  STO 13   35  RCL 05   36  REGSWAP   37  REGMOVE   38  XEQ IND 00   39  LBL 02   40  RCL IND 07   41  RCL 02   42  *   43  STO IND 08   44  3   45  /   46  ST+ IND 06   47  DSE 08   48  DSE 07   49  DSE 06   50  GTO 02   51  RCL 01 52  ST+ 06   53  ST+ 07   54  ST+ 08   55  RCL 02   56  LASTX   57  /   58  ST+ 20   59  XEQ IND 00   60  RCL 05   61  REGMOVE   62  LBL 03   63  RCL IND 07   64  RCL 02   65  *   66  STO IND 09   67  1.5   68  /   69  ST+ IND 06   70  DSE 09   71  DSE 07   72  DSE 06   73  GTO 03   74  RCL 01   75  ST+ 06   76  ST+ 07   77  ST+ 09   78  RCL 02   79  3   80  /   81  ST+ 20   82  XEQ IND 00   83  RCL 05   84  REGMOVE   85  LBL 04   86  RCL IND 07   87  RCL 02   88  *   89  STO IND 10   90  RCL IND 09   91  4   92  *   93  -   94  RCL IND 08   95  -   96  12   97  /   98  ST- IND 06   99  DSE 10 100  DSE 09 101  DSE 08 102  DSE 07 103  DSE 06 104  GTO 04 105  RCL 01 106  ST+ 06 107  ST+ 07 108  ST+ 08 109  ST+ 09 110  ST+ 10 111  RCL 02 112  3 113  / 114  ST- 20 115  XEQ IND 00 116  RCL 05 117  REGMOVE 118  LBL 05 119  RCL IND 07 120  RCL 02 121  * 122  STO IND 11 123  18 124  * 125  RCL IND 10 126  7 127  * 128  + 129  RCL IND 09 130  22 131  * 132  - 133  RCL IND 08 134  5 135  * 136  + 137  9.6 138  / 139  ST+ IND 06 140  DSE 11 141  DSE 10 142  DSE 09 143  DSE 08 144  DSE 07 145  DSE 06 146  GTO 05 147  RCL 01 148  ST+ 06 149  ST+ 07 150  ST+ 08 151  ST+ 09 152  ST+ 10 153  ST+ 11 154  RCL 02 155  2 156  / 157  ST+ 20 158  XEQ IND 00 159  RCL 05 160  REGMOVE 161  LBL 06 162  RCL IND 07 163  RCL 02 164  * 165  STO IND 12 166  5 167  / 168  RCL IND 11 169  + 170  RCL IND 10 171  4 172  / 173  - 174  RCL IND 08 175  18 176  * 177  RCL IND 09 178  55 179  * 180  - 181  60 182  / 183  + 184  2 185  / 186  ST+ IND 06 187  DSE 12 188  DSE 11 189  DSE 10 190  DSE 09 191  DSE 08 192  DSE 07 193  DSE 06 194  GTO 06 195  RCL 01 196  ST+ 06 197  ST+ 07 198  ST+ 08 199  ST+ 09 200  ST+ 10 201  ST+ 11 202  ST+ 12 203  RCL 02 204  6 205  / 206  RCL 13 207  + 208  STO 20 209  XEQ IND 00 210  RCL 05 211  REGMOVE 212  LBL 07 213  RCL IND 07 214  RCL 02 215  * 216  80 217  X<>Y 218  ST* Y 219  X<> IND 09 220  99 221  * 222  + 223  RCL IND 11 224  118 225  * 226  - 227  20 228  * 229  RCL IND 12 230  ST+ IND 09 231  128 232  * 233  + 234  RCL IND 10 235  ST+ IND 11 236  215 237  * 238  + 239  RCL IND 08 240  783 241  * 242  - 243  780 244  / 245  ST+ IND 06 246  DSE 12 247  DSE 11 248  DSE 10 249  DSE 09 250  DSE 08 251  DSE 07 252  DSE 06 253  GTO 07 254  RCL 01 255  ST+ 06 256  ST+ 07 257  ST+ 08 258  ST+ 09 259  ST+ 10 260  ST+ 11 261  ST+ 12 262  RCL 02 263  RCL 13 264  + 265  STO 20 266  XEQ IND 00 267  RCL 05 268  REGMOVE 269  LBL 08 270  RCL IND 07 271  RCL 02 272  * 273  RCL IND 08 274  + 275  13 276  * 277  RCL IND 09 278  32 279  * 280  + 281  RCL IND 11 282  55 283  * 284  + 285  .5 286  % 287  ST+ IND 06 288  DSE 11 289  DSE 09 290  DSE 08 291  DSE 07 292  DSE 06 293  GTO 08 294  RCL 01 295  ST+ 06 296  ST+ 07 297  ST+ 08 298  ST+ 09 299  ST+ 11 300  DSE 04 301  GTO 01 302  RCL 23 303  RCL 22 304  RCL 21 305  RCL 20 306  END

( 498 bytes / SIZE 21+8n )

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

Example:

dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1   ( whence N = 10 )
dz/dx = x.( y + z - u )         z(0) = 1
du/dx = x.y - z.u                u(0) = 2

LBL "T"      *                 RCL 21      -                   RCL 20      RCL 23       RTN            "T"  ASTO 00
RCL 21      *                 RCL 22      RCL 20        RCL 21      *
RCL 22      CHS           +                 *                  *                 -
RCL 23      STO 24      RCL 23       STO 25        RCL 22      STO 26

Initial values:     0   STO 20

1   STO 21    STO 22
2   STO 23

n = 3      STO 01    ( 3 functions )
h = 0.1   STO 02
N = 10   STO 03    XEQ "RK6"  >>>>      x   = 1                     = R20                      ( in 6mn28s )
RDN    y(1) = 0.258207889  = R21                                                                         y(1) = 0.258207906
RDN    z(1) = 1.157623948  = R22         the exact results, rounded to 9D are          z(1) = 1.157623981
RDN    u(1) = 0.842178329  = R23                                                                         u(1) = 0.842178312

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if needed )

c) A Runge-Kutta Method of Order 8

Data Registers:           •  R00 = subroutine name     ( Registers R00 to R03 , R14 to R48 and R50 to R50+n are to be initialized before executing "RK8" )

•  R01 = n  = number of equations = number of functions        R04 thru R13: temp        R49: temp
•  R02 = h  = stepsize
•  R03 = N = number of steps                                            •  R14 thru R48 = Constants

•  R50 = x0
•  R51 = y1(x0)                                                                           R51+n thru R50+9n: temp
•  R52 =  y2(x0)
...............         ( initial values )

•  R50+n = yn(x0)

•  R14 =  0.8961811934                 •  R23 =  0.1996369936                    •  R32 =  0.9576053440                    •  R41 =  0.9401094520
•  R15 = -0.2117115009                •  R24 =  0.09790045616                  •  R33 =  -0.6325461607                   •  R42 = -2.229158210
•  R16 =  0.8273268354                 •  R25 =  0.3285172131                    •  R34 =  0.1527777778                    •  R43 =  7.553840442
•  R17 =  0.06514850915               •  R26 = -0.3497052863                    •  R35 = -0.009086961101               •  R44 = -7.164951553
•  R18 =  0.5766714727                 •  R27 = -0.03302551131                  •  R36 =  1.062498447                      •  R45 =  2.451380432
•  R19 =  0.1855068535                 •  R28 =  0.1289862930                     •  R37 = -1.810863083                     •  R46 =  0.05
•  R20 =  0.3865246267                 •  R29 =  0.1726731646                     •  R38 =  2.031083139                     •  R47 =   0.2722222222
•  R21 = -0.4634553896                •  R30 =  -0.01186868389                  •  R39 = -0.6379313502                   •  R48 =  0.3555555556
•  R22 =  0.3772937693                 •  R31 =  0.002002165993                 •  R40 =  -0.5512205631

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

-Line 459 is a three-byte GTO 01

 01  LBL "RK8"   02  51.051   03  RCL 01   04   E6   05  /   06  +   07  RCL 01   08  8   09  *   10  +   11  STO 05   12  RCL 01   13  RCL 01   14  RCL 01   15  50.05   16  +   17  STO 06   18  +   19  STO 07   20  +   21  STO 08   22  +   23  STO 09   24  +   25  STO 10   26  +   27  STO 11   28  +   29  STO 12   30  +   31  STO 13   32  RCL 03   33  STO 49   34  GTO 01   35  LBL 00   36  XEQ IND 00   37  RCL 05   38  REGMOVE   39  RTN   40  LBL 01   41  RCL 50   42  STO 04   43  RCL 05   44  REGSWAP   45  REGMOVE   46  XEQ IND 00   47  LBL 02   48  RCL IND 07   49  RCL 02   50  *   51  STO IND 08   52  2   53  /   54  ST+ IND 06   55  DSE 08   56  DSE 07   57  DSE 06   58  GTO 02   59  RCL 01   60  ST+ 06   61  ST+ 07   62  ST+ 08   63  RCL 02   64  LASTX   65  /   66  ST+ 50   67  XEQ 00   68  LBL 03   69  RCL IND 07   70  RCL 02   71  *   72  STO IND 09   73  RCL IND 08   74  +   75  4   76  /   77  ST+ IND 06   78  DSE 09 79  DSE 08   80  DSE 07   81  DSE 06   82  GTO 03   83  RCL 01   84  ST+ 06   85  ST+ 07   86  ST+ 08   87  ST+ 09   88  XEQ 00   89  LBL 04   90  RCL IND 07   91  RCL 02   92  *   93  STO IND 10   94  RCL 14   95  *   96  RCL IND 09   97  RCL 15   98  *   99  + 100  RCL IND 08 101  7 102  / 103  + 104  ST+ IND 06 105  DSE 10 106  DSE 09 107  DSE 08 108  DSE 07 109  DSE 06 110  GTO 04 111  RCL 01 112  ST+ 06 113  ST+ 07 114  ST+ 08 115  ST+ 09 116  ST+ 10 117  RCL 02 118  RCL 16 119  * 120  RCL 04 121  + 122  STO 50 123  XEQ 00 124  LBL 05 125  RCL IND 07 126  RCL 02 127  * 128  STO IND 11 129  RCL 17 130  * 131  RCL IND 10 132  RCL 18 133  * 134  + 135  RCL IND 08 136  RCL 19 137  * 138  + 139  ST+ IND 06 140  DSE 11 141  DSE 10 142  DSE 08 143  DSE 07 144  DSE 06 145  GTO 05 146  RCL 01 147  ST+ 06 148  ST+ 07 149  ST+ 08 150  ST+ 10 151  ST+ 11 152  XEQ 00 153  LBL 06 154  RCL IND 07 155  RCL 02 156  * 157  STO IND 12 158  RCL 20 159  * 160  RCL IND 11 161  RCL 21 162  * 163  + 164  RCL IND 10 165  RCL 22 166  * 167  + 168  RCL IND 08 169  RCL 23 170  * 171  + 172  ST+ IND 06 173  DSE 12 174  DSE 11 175  DSE 10 176  DSE 08 177  DSE 07 178  DSE 06 179  GTO 06 180  RCL 01 181  ST+ 06 182  ST+ 07 183  ST+ 08 184  ST+ 10 185  ST+ 11 186  ST+ 12 187  RCL 02 188  2 189  / 190  RCL 04 191  + 192  STO 50 193  XEQ 00 194  LBL 07 195  RCL IND 07 196  RCL 02 197  * 198  STO IND 13 199  RCL 24 200  * 201  RCL IND 12 202  RCL 25 203  * 204  + 205  RCL IND 11 206  RCL 26 207  * 208  + 209  RCL IND 10 210  RCL 27 211  * 212  + 213  RCL IND 08 214  RCL 28 215  * 216  + 217  ST+ IND 06 218  DSE 13 219  DSE 12 220  DSE 11 221  DSE 10 222  DSE 08 223  DSE 07 224  DSE 06 225  GTO 07 226  RCL 01 227  ST+ 06 228  ST+ 07 229  ST+ 08 230  ST+ 10 231  ST+ 11 232  ST+ 12 233  ST+ 13 234  RCL 02 235  RCL 29 236  * 237  RCL 04 238  + 239  STO 50 240  XEQ 00 241  LBL 08 242  RCL IND 07 243  RCL 02 244  * 245  STO IND 09 246  9 247  / 248  RCL IND 13 249  RCL 30 250  * 251  + 252  RCL IND 12 253  RCL 31 254  * 255  + 256  RCL IND 08 257  14 258  / 259  + 260  ST+ IND 06 261  DSE 13 262  DSE 12 263  DSE 09 264  DSE 08 265  DSE 07 266  DSE 06 267  GTO 08 268  RCL 01 269  ST+ 06 270  ST+ 07 271  ST+ 08 272  ST+ 09 273  ST+ 12 274  ST+ 13 275  XEQ 00 276  LBL 09 277  RCL IND 07 278  RCL 02 279  * 280  STO IND 10 281  RCL 32 282  * 283  RCL IND 09 284  RCL 33 285  * 286  + 287  RCL IND 13 288  RCL 34 289  * 290  + 291  RCL IND 12 292  RCL 35 293  * 294  + 295  RCL IND 08 296  32 297  / 298  + 299  ST+ IND 06 300  DSE 13 301  DSE 12 302  DSE 10 303  DSE 09 304  DSE 08 305  DSE 07 306  DSE 06 307  GTO 09 308  RCL 01 309  ST+ 06 310  ST+ 07 311  ST+ 08 312  ST+ 09 313  ST+ 10 314  ST+ 12 315  ST+ 13 316  RCL 02 317  2 318  / 319  RCL 04 320  + 321  STO 50 322  XEQ 00 323  LBL 10 324  RCL IND 07 325  RCL 02 326  * 327  STO IND 11 328  RCL 36 329  * 330  RCL IND 10 331  RCL 37 332  * 333  + 334  RCL IND 09 335  RCL 38 336  * 337  + 338  RCL IND 13 339  RCL 39 340  * 341  + 342  RCL IND 12 343  9 344  / 345  + 346  RCL IND 08 347  14 348  / 349  + 350  ST+ IND 06 351  DSE 13 352  DSE 12 353  DSE 11 354  DSE 10 355  DSE 09 356  DSE 08 357  DSE 07 358  DSE 06 359  GTO 10 360  RCL 01 361  ST+ 06 362  ST+ 07 363  ST+ 08 364  ST+ 09 365  ST+ 10 366  ST+ 11 367  ST+ 12 368  ST+ 13 369  RCL 02 370  RCL 16 371  * 372  RCL 04 373  + 374  STO 50 375  XEQ 00 376  LBL 11 377  RCL IND 07 378  RCL 02 379  * 380  X<> IND 12 381  RCL 40 382  * 383  RCL IND 12 384  RCL 41 385  * 386  + 387  RCL IND 11 388  RCL 42 389  * 390  + 391  RCL IND 10 392  RCL 43 393  * 394  + 395  RCL IND 09 396  RCL 44 397  * 398  + 399  RCL IND 13 400  RCL 45 401  * 402  + 403  ST+ IND 06 404  DSE 13 405  DSE 12 406  DSE 11 407  DSE 10 408  DSE 09 409  DSE 07 410  DSE 06 411  GTO 11 412  RCL 01 413  ST+ 06 414  ST+ 07 415  ST+ 09 416  ST+ 10 417  ST+ 11 418  ST+ 12 419  ST+ 13 420  RCL 02 421  RCL 04 422  + 423  STO 50 424  XEQ 00 425  LBL 12 426  RCL IND 07 427  RCL 02 428  * 429  RCL IND 08 430  + 431  RCL 46 432  * 433  RCL IND 10 434  RCL IND 12 435  + 436  RCL 47 437  * 438  + 439  RCL IND 11 440  RCL 48 441  * 442  + 443  ST+ IND 06 444  DSE 12 445  DSE 11 446  DSE 10 447  DSE 08 448  DSE 07 449  DSE 06 450  GTO 12 451  RCL 01 452  ST+ 06 453  ST+ 07 454  ST+ 08 455  ST+ 10 456  ST+ 11 457  ST+ 12 458  DSE 49 459  GTO 01  460  RCL 53 461  RCL 52 462  RCL 51 463  RCL 50 464  END

( 765 bytes / SIZE 51+9n )

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

Example:

dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1  and  N = 10
dz/dx = x.( y + z - u )         z(0) = 1
du/dx = x.y - z.u                u(0) = 2

LBL "T"      *                 RCL 51      -                   LASTX      RCL 53       RTN            "T"  ASTO 00
RCL 51      *                 RCL 52      RCL 50        RCL 51      *
RCL 52      CHS           +                 *                  *                 -
RCL 53      STO 54      RCL 53       STO 55        RCL 52      STO 56

Initial values:     0   STO 50

1   STO 51    STO 52
2   STO 53

n = 3      STO 01    ( 3 functions )
h = 0.1   STO 02
N = 10   STO 03    XEQ "RK8"  >>>>      x   = 1                     = R50                      ( in 9mn31s )
RDN    y(1) = 0.258207906  = R51                                                                         y(1) = 0.258207906
RDN    z(1) = 1.157623979  = R52         the exact results, rounded to 9D are          z(1) = 1.157623981
RDN    u(1) = 0.842178313  = R53                                                                         u(1) = 0.842178312

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if needed )
-If you replace lines 431 thru 442 by  9  *  RCL IND 10  RCL IND 12  +  49  *  +  RCL IND 11  64  *  +  PI  R-D  /
registers  R46-R47-R48  will be unused.

d) Runge-Kutta-Fehlberg 4-5 Method

-The following program uses a 4th-order formula to compute the solution,
and the difference between this 4th-order formula and a 5th-order formula provides an error estimate.

Data Registers:           •  R00 = subroutine name       ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RKF" )

•  R01 = n  = number of equations = number of functions         R04 thru R14: temp        R15 thru R19: unused
•  R02 = h  = stepsize
•  R03 = N = number of steps

•  R20 = x0                                         at the end:
•  R21 = y1(x0)                                     R21+n = err(y1)                                      R21+n thru R20+8n: temp
•  R22 =  y2(x0)                                    R22+n = err(y2)
...............                                      ......................

•  R20+n = yn(x0)                                R20+2n = err(yn)

Where  err(yi) is an error-estimate of yi    ( during the calculations, these error-estimates are moved in registers R21+2n .... R20+3n )

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

-If you don't have an HP-41 CX, replace line 31 by   0   LBL 00   STO IND Y   ISG Y   GTO 00
-Lines 290 & 298 are three-byte GTOs

 01  LBL "RKF"    02  21.021   03  RCL 01   04   E6   05  /   06  +   07  RCL 01   08  7   09  *   10  +   11  STO 05   12  FRC   13  RCL 01   14  .1   15  %   16  +   17  +   18  RCL 01   19  +   20  21   21  +   22  STO 06   23  RCL 01   24  ST+ X   25   E3   26  /   27  RCL 01   28  +   29  21.02   30  +   31  CLRGX   32  RCL 01   33  RCL 01   34  LASTX   35  1   36  -   37  +   38  STO 07   39  +   40  STO 08   41  +   42  STO 13   43  +   44  STO 09   45  +   46  STO 10   47  +   48  STO 11   49  +   50  STO 12 51  LBL 10   52  RCL 03   53  STO 04   54  RCL 06   55  REGSWAP   56  LBL 01   57  RCL 20   58  STO 14   59  RCL 05   60  REGSWAP   61  REGMOVE   62  XEQ IND 00   63  LBL 02   64  RCL IND 08   65  RCL 02   66  *   67  STO IND 09   68  ST+ X   69  9   70  /   71  ST+ IND 07   72  DSE 09   73  DSE 08   74  DSE 07   75  GTO 02   76  RCL 01   77  ST+ 07   78  ST+ 08   79  ST+ 09   80  RCL 02   81  ST+ X   82  LASTX   83  /   84  ST+ 20   85  XEQ IND 00   86  RCL 05   87  REGMOVE   88  LBL 03   89  RCL IND 08   90  RCL 02   91  *   92  4   93  /   94  STO IND 10   95  RCL IND 09   96  12   97  /   98  +   99  ST+ IND 07 100  DSE 10 101  DSE 09 102  DSE 08 103  DSE 07 104  GTO 03 105  RCL 01 106  ST+ 07 107  ST+ 08 108  ST+ 09 109  ST+ 10 110  RCL 02 111  9 112  / 113  ST+ 20 114  XEQ IND 00 115  RCL 05 116  REGMOVE 117  LBL 04 118  RCL IND 08 119  RCL 02 120  * 121  STO IND 11 122  270 123  * 124  RCL IND 10 125  972 126  * 127  - 128  RCL IND 09 129  69 130  * 131  + 132  128 133  / 134  ST+ IND 07 135  DSE 11 136  DSE 10 137  DSE 09 138  DSE 08 139  DSE 07 140  GTO 04 141  RCL 01 142  ST+ 07 143  ST+ 08 144  ST+ 09 145  ST+ 10 146  ST+ 11 147  RCL 02 148  2.4 149  / 150  ST+ 20 151  XEQ IND 00 152  RCL 05 153  REGMOVE 154  LBL 05 155  RCL IND 08 156  RCL 02 157  * 158  64 159  * 160  STO IND 12 161  RCL IND 09 162  85 163  * 164  - 165  60 166  / 167  RCL IND 10 168  RCL IND 11 169  5 170  / 171  - 172  27 173  * 174  + 175  ST+ IND 07 176  DSE 12 177  DSE 11 178  DSE 10 179  DSE 09 180  DSE 08 181  DSE 07 182  GTO 05 183  RCL 01 184  ST+ 07 185  ST+ 08 186  ST+ 09 187  ST+ 10 188  ST+ 11 189  ST+ 12 190  RCL 02 191  ST+ 14 192  RCL 14 193  STO 20 194  XEQ IND 00 195  RCL 05 196  REGMOVE 197  LBL 06 198  RCL IND 08 199  RCL 02 200  * 201  15 202  * 203  ST+ IND 12 204  RCL IND 12 205  RCL IND 11 206  351 207  * 208  + 209  RCL IND 10 210  540 211  * 212  - 213  RCL IND 09 214  65 215  * 216  + 217  432 218  / 219  ST+ IND 07 220  DSE 12 221  DSE 11 222  DSE 10 223  DSE 09 224  DSE 08 225  DSE 07 226  GTO 06 227  RCL 01 228  ST+ 07 229  ST+ 08 230  ST+ 09 231  ST+ 10 232  ST+ 11 233  ST+ 12 234  RCL 02 235  6 236  / 237  ST- 20 238  XEQ IND 00 239  RCL 05 240  REGMOVE 241  LBL 07 242  RCL IND 08 243  RCL 02 244  * 245  72 246  * 247  RCL IND 12 248  - 249  RCL IND 11 250  9 251  * 252  + 253  RCL IND 09 254  ST+ X 255  - 256  3 257  1/X 258  % 259  ST- IND 13 260  RCL IND 12 261  RCL IND 11 262  81 263  * 264  + 265  PI 266  R-D 267  / 268  RCL IND 09 269  9 270  / 271  + 272  ST+ IND 07 273  DSE 13 274  DSE 12 275  DSE 11 276  DSE 09 277  DSE 08 278  DSE 07 279  GTO 07 280  RCL 01 281  ST+ 07 282  ST+ 08 283  ST+ 09 284  ST+ 11 285  ST+ 12 286  ST+ 13 287  RCL 14 288  STO 20 289  DSE 04 290  GTO 01 291  RCL 06 292  REGSWAP 293  RCL 23 294  RCL 22 295  RCL 21 296  RCL 20 297  RTN 298  GTO 10 299  END

( 481 bytes / SIZE 21+8n )

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

Example:

dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1  &  N = 10
dz/dx = x.( y + z - u )         z(0) = 1
du/dx = x.y - z.u                u(0) = 2

LBL "T"      *                 RCL 21      -                   RCL 20      RCL 23       RTN            "T"  ASTO 00
RCL 21      *                 RCL 22      RCL 20        RCL 21      *
RCL 22      CHS           +                 *                  *                 -
RCL 23      STO 24      RCL 23       STO 25        RCL 22      STO 26

Initial values:     0   STO 20

1   STO 21    STO 22
2   STO 23

n = 3      STO 01    ( 3 functions )
h = 0.1   STO 02
N = 10   STO 03    XEQ "RKF"  >>>>      x   = 1                     = R20                      ( in 5mn40s )
RDN    y(1) = 0.258207319  = R21                                                              err(y) = -603 E-9 = R24
RDN    z(1) = 1.157624973  = R22                and the error estimates            err(z) = 1062 E-9 = R25
RDN    u(1) = 0.842178529  = R23                                                              err(u) = 1768 E-9 = R26

-Actually, the true errors are    -587 E-9  for y(1)     992 E-9  for z(1)   217 E-9  for u(1)

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if needed )

-If you want to use the 5th-order formula to compute yi(x) , replace line 259 with ST+ IND 07   ST+ IND 13  but the errors will be probably overestimated.
-In the above example, it gives:

y(1) = 0.258207897
z(1) = 1.157624056
u(1) = 0.842178340

e)  Bulirsch-Stoer Method

-This program uses the formulae given in "The Bulirsch-Stoer Method for the HP-41" and - first of all - in "Numerical Recipes"

Data Registers:           •  R00 = subroutine name       ( Registers R00-R01-R02 and R20 thru R20+n are to be initialized before executing "BST" )

•  R01 = n  = number of equations = number of functions         R03 thru R18: temp        R19: unused
•  R02 = tolerance = a small positive number
that specifies the desired accuracy

•  R20 = x0
•  R21 = y1(x0)                                                                           R21+n thru R20+14n: temp
•  R22 =  y2(x0)
...............         ( initial values )

•  R20+n = yn(x0)
Flags:  F09-F10
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+n, ... , R20+2n  respectively
with  x, y1, ... , yn in regiters R20, R21, ... , R20+n

-Lines 194-245-252-262-266  are three-byte  GTOs

 01  LBL "BST"   02  STO 03   03  1   04  STO 04   05  9   06  STO 05   07  21.021   08  RCL 01   09   E6   10  /   11  +   12  RCL 01   13  5   14  *   15  +   16  STO 08   17  RCL 01   18  ST+ X   19  ST- Y   20  .1   21  %   22  -   23  -   24  STO 16          25  20.02   26  RCL 01   27  +   28  STO 10   29  INT   30  .1   31  %   32  RCL 01   33  STO T   34  +   35  +   36  STO 11   37  +   38  STO 12   39  +   40  STO 13   41  +   42  STO 14   43  GTO 01   44  LBL 00   45  2 46  ST/ 04   47  GTO 02   48  LBL 01   49  RCL 20   50  STO 15   51  RCL 08   52  REGSWAP   53  REGMOVE   54  XEQ IND 00   55  RCL 16   56  REGMOVE   57  6   58  RCL 05   59  X>Y?   60  GTO 02   61  RCL 04   62  ST+ 04   63  LBL 02   64  SF 09   65  RCL 04   66  STO 17   67  ABS   68  RCL 03   69  RCL 15   70  -   71  ABS   72  X<=Y?   73  CF 09   74  X>Y?   75  X<>Y   76  LASTX   77  SIGN   78  *   79  STO 04   80  SF 10   81  CLX   82  STO 05   83  RCL 08   84  INT   85  STO 09   86  LBL 03   87  RCL 08   88  REGMOVE    89  2   90  ST+ 05 91  18   92  RCL 05   93  X=Y?   94  GTO 00   95  RCL 01   96  ST+ 09   97  SIGN   98  -   99   E3 100  / 101  ISG X 102  STO 07 103  RCL 04 104  RCL 05 105  / 106  STO 06 107  LBL 04 108  RCL IND 12 109  RCL 06 110  * 111  RCL IND 10 112  STO IND 13 113  + 114  STO IND 14 115  DSE 14 116  DSE 13 117  DSE 12 118  DSE 10 119  GTO 04 120  RCL 01 121  ST+ 10 122  ST+ 12 123  ST+ 13 124  ST+ 14 125  LBL 05 126  RCL 08 127  RCL 01 128  - 129  REGMOVE 130  RCL 06 131  RCL07 132  INT 133  * 134  RCL 15 135  + 136  STO 20 137  XEQ IND 00 138  LBL 06 139  RCL IND 11 140  RCL 06 141  ST+ X 142  * 143  RCL IND 14 144  X<> IND 13 145  + 146  STO IND 14 147  DSE 14 148  DSE 13 149  DSE 11 150  GTO 06 151  RCL 01 152  ST+ 11 153  ST+ 13 154  ST+ 14 155  ISG 07 156  GTO 05 157  RCL 08 158  RCL 01 159  - 160  REGMOVE 161  RCL 04 162  RCL 15 163  + 164  STO 20 165  XEQ IND 00 166  LBL 07 167  RCL IND 11 168  RCL 06 169  * 170  RCL IND 13 171  + 172  ST+ IND 10 173  2 174  ST/ IND 10 175  SIGN 176  ST- 11 177  ST- 13 178  DSE 10 179  GTO 07 180  21 181  RCL 09 182   E3 183  / 184  + 185  RCL 01 186  ST+ 10 187  ST+ 11 188  ST+ 13 189   E6 190  / 191  + 192  REGMOVE 193  FS?C 10 194  GTO 03   195  RCL 09 196  RCL 01 197  + 198  DSE X 199   E3 200  / 201  RCL 09 202  + 203  STO 06 204  CLX 205  STO 18 206  LBL 08 207  RCL 05 208  2 209  ST- Y 210  / 211  STO 07 212  LBL 09 213  RCL 06 214  RCL 07 215  RCL 01 216  * 217  - 218  RCL IND 06 219  ENTER^ 220  X<> IND Z 221  STO Z 222  - 223  RCL 05 224  X^2 225  ST* Y 226  RCL 07 227  ST+ X 228  X^2 229  - 230  / 231  + 232  STO IND 06 233  DSE 07 234  GTO 09 235  LASTX 236  ABS 237  RCL 18 238  XY 240  STO 18 241  ISG 06 242  GTO 08 243  RCL 02 244  X

( 395 bytes / SIZE 21+14n )

 STACK INPUTS OUTPUTS T / y3(x1) Z / y2(x1) Y / y1(x1) X x1 x1

Example:

dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  a tolerance = 10 -7
dz/dx = x.( y + z - u )         z(0) = 1
du/dx = x.y - z.u                u(0) = 2

LBL "T"      *                 RCL 21      -                   LASTX      RCL 23       RTN            "T"  ASTO 00
RCL 21      *                 RCL 22      RCL 20        RCL 21      *
RCL 22      CHS           +                 *                  *                 -
RCL 23      STO 24      RCL 23       STO 25        RCL 22      STO 26

Initial values:     0   STO 20

1   STO 21    STO 22
2   STO 23

n  =  3     STO 01    ( 3 functions )
tol = E-7  STO 02

1   XEQ "BST"  >>>>    x   =  1                    = R20                      ( in 12mn06s )
RDN    y(1) = 0.258207909  = R21                                                                         y(1) = 0.258207906
RDN    z(1) = 1.157623986  = R22         the exact results, rounded to 9D are          z(1) = 1.157623981
RDN    u(1) = 0.842178304  = R23                                                                         u(1) = 0.842178312

-The long execution time is due to the fact that the method is first used with H = 1
but the desired accuracy is not satisfied. So the stepsize is divided by 2.

-To compute  y(2)  z(2)  u(2)  key in  2  R/S   it yields              ( in 6m46s )

y(2) = 0.106363294                                                 y(2) = 0.106363329
z(2) = 3.886706181         the exact results are          z(2) = 3.886706159         ( rounded to 9D )
u(2) = 0.196515847                                                 u(2) = 0.196515847

Notes:

-The initial stepsize is always +/-1 ( line 03 )
-If you want to choose your own initial stepsize, place it in Y-register after replacing line 03 by X<>Y
-The modified midpoint rule ( and Richarson extrapolation ) are used with  h/2 , h/4 , h/6 , ..... , h/16 ( at most )  until the desired accuracy is satisfied.
-If this is not satisfied, even with h/16, the stepsize is divided by 2 ( LBL 00 )
-If - for instance - you want to stop the extrapolation with  h/14, replace line 91 by 16 ( instead of 18 )
-Thus, the SIZE may be reduced to 21+13n
-On the other hand, if you want to continue the extrapolation with  h/18, replace line 91 by 20
-In this case, the SIZE must be at least 21+15n

2°) Systems of second-order Conservative Equations

-The 3 following programs use the 4th-order Runge-Kutta formula:

y(x+h) = y(x) + h ( y'(x) + k1/6 + 2k2/3 )
y'(x+h) = y'(x) + k1/6 + 2k2/3 + k3/6                 where  k1 = h.f (x,y)  ;  k2 = h.f(x+h/2,y+h.y'/2+h.k1/8)  ;  k3 = h.f(x+h,y+h.y'+h.k2/2)

a)  Systems of 2 Equations

-"RKS2" 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:  subroutine name                           ( Registers R00 thru R07  are to be initialized before executing "RKS2" )

•   R01 = x0      •   R04 = h/2 = half of the stepsize       •   R06 = y'0          R08 to R12: temp
•   R02 = y0      •   R05 = N = number of steps            •   R07 = z'0
•   R03 = z0
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.

-Line 96 is a three-byte  GTO 01

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

( 139 bytes / SIZE 013 )

 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

LBL "T"      CHS        *                     "T" ASTO 00
RDN          ST* Z       RTN
STO Z        -
X<>Y        R^

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

XEQ "RKS2"  >>>>     x   =  1                                                                                                  ( in 41 seconds )
RDN   y(1) = 1.531358015     and   R06 = y'(1) = -2.312838895
RDN   z(1) = 2.620254480              R07 = z'(1) =  2.941751649

-With  h = 0.05  it yields   y(1) =  1.531356736      y'(1) = -2.312840085
z(1) =  2.620254295      z'(1) =  2.941748608

-Actually, the exact results - rounded to 9D - are

y(1) =  1.531356646      y'(1) = -2.312840137
z(1) =  2.620254281      z'(1) =  2.941748399

-To continue the calculations, simply press R/S  ( after changing h/2 & N in registers R04 & R05 if needed )

b)  Systems of 3 Equations

-"RKS3" solves the system:

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

Data Registers:    •   R00:  subroutine name                           ( Registers R00 thru R09  are to be initialized before executing "RKS3" )

•   R01 = x0      •   R05 = h/2 = half of the stepsize       •   R07 = y'0
•   R02 = y0      •   R06 = N = number of steps            •   R08 = z'0
•   R03 = z0                   R10 to R16: temp                  •   R09 = u'0
•   R04 = u0
Flags: /
Subroutine:  a program which calculates f(x;y;z;u) in Z-register , g(x;y;z;u) in Y-register and  h(x;y;z;u) in X-register
assuming x is in X-register , y is in Y-register , z is in Z-register and u is in T-register upon entry.

-Line 136 is a three-byte  GTO 01

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

( 194 bytes / SIZE 017 )

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

Example:

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

LBL "T"       R^               RCL Z         X<> Z         LASTX                 "T"  ASTO 00
R^               ST* Y         R^                ST+ Y         *
CHS            ST+ L         ST+ L          ST* Z          X<>Y
STO L         CLX           ST* Y          X<> T         RTN

0     STO 01
1     STO 02  STO 03  STO 07  STO 08  STO 09
2     STO 04
0.05  STO 05
10    STO 06    XEQ "RKS3"  >>>>     x   = 1                                                                                             ( in 65 seconds )
RDN   y(1) = 0.439528419           y'(1) = -2.101120400 = R07
RDN   z(1) = 2.070938499           z'(1) =   1.269599239 = R08
RDN   u(1) = 1.744522976           u'(1) = -1.704232092 = R09

-With  h = 0.05  it yields

y(1) = 0.439524393          y'(1) = -2.101122784
z(1) = 2.070940521          z'(1) =   1.269597110
u(1) = 1.744524843          u'(1) = -1.704234567

-Actually, the exact results rounded to 9D are:

y(1) = 0.439524100          y'(1) = -2.101122880
z(1) = 2.070940654          z'(1) =   1.269596950
u(1) = 1.744524964          u'(1) = -1.704234756

-Press R/S  to continue the calculations ( after changing h & N in registers R05 & R06 if needed )
-The subroutine may be difficult to write, but "RKS3" is faster than "RKS" below, for 3 equations.

c)  Systems of n Equations

Data Registers:           •  R00 = subroutine name       ( Registers R00 thru R03 and R10 thru R10+2n are to be initialized before executing "RKS" )

•  R01 = n  = number of equations = number of functions         R04 thru R09: temp
•  R02 = h  = stepsize
•  R03 = N = number of steps

•  R10 = x0
•  R11 = y1(x0)                                  •  R11+n = y'1(x0)                                      R11+2n thru R10+6n: temp
•  R12 =  y2(x0)                                 •  R12+n = y'2(x0)
...............                                      ......................

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

( during the calculations, y'i  are moved in registers R11+2n .... R10+3n )

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

-Line 140 is a three-byte  GTO 01

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

( 225 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:

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

LBL "T"      *                 RCL 11      -                   RCL 10      RCL 13       RTN            "T"  ASTO 00
RCL 11      *                 RCL 12      RCL 10        RCL 11      *
RCL 12      CHS           +                 *                  *                 -
RCL 13      STO 14      RCL 13       STO 15        RCL 12      STO 16

Initial values:     0   STO 10

1   STO 11    STO 12    STO 14    STO 15    STO 16
2   STO 13

n = 3      STO 01    ( 3 functions )
h = 0.1   STO 02
N = 10   STO 03    XEQ "RKS"  >>>>    x   =  1                     = R10                                                                ( in 2mn19s )
RDN   y(1) = 0.439528419  = R11          y'(1) = -2.101120401 = R14
RDN   z(1) = 2.070938499  = R12          z'(1) =   1.269599239 = R15
RDN   u(1) = 1.744522977  = R13          u'(1) = -1.704232091 = R16

-Press R/S  to continue the calculations ( after changing h & N in registers R02 & R03 if needed )

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[3]  F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[4]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4