# hp41programs

Systems Linear and non-Linear Systems for the HP-41

Overview

1°) Linear Systems

a) Program#1
b) Program#2
c) Program#3
d) Underdetermined Systems

d-1) Program#1
d-2) Program#2

e) Overdetermined Systems

e-1) Program#1
e-2) Program#2

f) Tridiagonal Systems

2°) Non-Linear Systems

a) 1 Equation in 1 unknown

a-1)  Secant Method
a-3)  Ridders' Method
a-4)  Discrete Data ( 5 points )

b) 2 Equations in 2 unknowns
c) 3 Equations in 3 unknowns
d) N Equations in N unknowns

d-1) Program#1
d-2) Program#2

3°) A short In/Out Routine

** To solve non-linear systems of equations ( including with complex unknowns ).

1°) Linear Systems

a) Program#1

"LS" allows you to solve linear systems, including overdetermined and underdetermined systems.
You can also invert up to a 12x12 matrix.
It is essentially equivalent to the RREF function of the HP-48 ( a "little" slower, I grant you that! ):
Its objective is to reduce the matrix on the upper left to a diagonal form with only ones in the diagonal.
The determinant of this matrix is also computed and stored in register R00.
( if there are more rows than columns, R00 is not always equal to the determinant of the upper left matrix because of rows exchanges )

If flag F01 is clear, Gaussian elimination with partial pivoting is used.
If flag F01 is set, the pivots are the successive elements of the diagonal:
It's sometimes useful for matrices like Pascal's matrices of high order.
They are extremely troublesome and many roundoff errors can occur.
But if you set flag F01, all the coefficients will be computed with no roundoff error at all, because all the pivots will be ones!

One advantage of this program is that you can choose the beginning data register - except R00 - ( this feature is used to solve non-linear systems too ):
You store the first coefficient into Rbb , ... , up to the last one into Ree ( column by column )             ( bb > 00 )

Then you key in  bbb.eeerr  ENTER^
p          XQ "LS"    and the system will be solved.    ( bbb.eeerr ends up in L-register )

where r is the number of rows of the matrix
and  p  is a small non-negative number that you choose to determine tiny elements: if a pivot is smaller than p it will be considered to be zero.

Don't interrupt "LS" because registers P and Q are used ( there is no risk of crash, but their contents could be altered )
If you're not familiar to synthetic programming, you can replace status registers M N O P Q by R01 R02 R03 R04 R05
In that case, you must choose b > 5.

-Line 121 is a three-byte  GTO 01

 01  LBL "LS"    02  STO M   03  X<>Y   04  ENTER^   05  INT   06  LASTX   07  FRC   08  ISG X   09  INT   10  STO O   11  RCLY   12  +   13  DSE X   14   E3   15  /   16  +   17  X<> O          18  .1   19  %   20  +   21  STO Q   22  SIGN   23  STO 00   24  X<>Y   25  STO P    26  LBL 01 27  STO N   28  FS? 01   29  GTO 04   30  INT   31  RCL O   32  FRC   33  +   34  ENTER^   35  ENTER^   36  CLX   37  LBL 02         38  RCL IND Z   39  ABS   40  X>Y?   41  STO Z   42  X>Y?   43  +   44  RDN   45  ISG Z   46  GTO 02   47  RCL N   48  ENTER^   49  FRC   50  R^   51  INT   52  + 53  X=Y?   54  GTO 04   55  LBL 03   56  RCL IND X   57  X<> IND Z   58  STO IND Y   59  ISG Y   60  RDN   61  ISG Y   62  GTO 03   63  RCL 00   64  CHS   65  STO 00   66  LBL 04   67  RCL M   68  RCL IND N   69  ST* 00   70  ABS   71  X<=Y?   72  GTO 09   73  RCL N   74  LASTX   75  LBL 05   76  ST/ IND Y   77  ISG Y   78  GTO 05 79  LBL 06   80  RCL N   81  ENTER^   82  FRC   83  RCL O   84  INT   85  +   86  X=Y?   87  GTO 08   88  RCL IND X   89  SIGN   90  RDN   91  LBL 07   92  RCL IND Y   93  LASTX   94  *   95  ST- IND Y   96  ISG Y   97  RDN   98  ISG Y   99  GTO 07 100  LBL 08 101  ISG O 102  GTO 06 103  RCL Q 104  INT 105  ST- O 106  LBL 09        107  RCL Q 108  ST+ O 109  X^2 110  RCL P 111  INT 112  ENTER^ 113  SIGN 114  ST+ N 115  - 116  + 117  RCL N  118  X>Y? 119  CLX 120  ISG X 121  GTO 01 122  X<> P 123  SIGN 124  RCL 00 125  CLA 126  END

( 189 bytes / SIZE eee+1 )

 STACK INPUTS OUTPUTS Y bbb.eeerr 1 X p determinant L / bbb.eeerr

Example1:     Find the solution of the system:

2x + 3y + 5z + 4t = 39
-4x + 2y + z + 3t  = 15
3x  -  y + 2z + 3t = 19
5x + 7y - 3z + 2t = 18

If you choose to store the first element into R01, you have to store these 20 numbers:

2  3  5  4  39                  R01  R05  R09  R13  R17
-4  2  1  3  15        in        R02  R06  R10  R14  R18      respectively           ( you can use "IO" at the bottom of this page )
3 -1  2  3  19                  R03  R07  R11  R15  R19
5  7 -3  2  18                  R04  R08  R12  R16  R20

The first register is R01, the last register is R20 and there are 4 rows, therefore the control number of the matrix is  1.02004
If you choose p = 10-7  key in:

CF 01
1.02004 ENTER^
E-7    XEQ "LS"   and  you obtain  840 in X-register and in R00:  it is the determinant of the 4x4 matrix on the left.              ( execution time = 31s )

Registers R01 thru R16 now contains the unit matrix and registers R17 thru R20 contains the solution  x = 1 ; y = 2 ; z = 3 ; t = 4
Thus, the array has been changed into:

1  0  0  0  1
0  1  0 0  2              the solution is the last column
0  0  1  0  3              of the new matrix.
0  0  0  1  4

Example2:    Solve the system:

5x + y + z  = 8
4x - y + 2z = 13
x + 2y - z  = -5
7x - 4y + 5z = 31

-Suppose we need to store the first element into  R11 , then we store these 16 numbers

5   1  1   8                       R11  R15  R19  R23
4  -1  2  13         in          R12  R16  R20  R24          respectively
1   2  -1  -5                     R13  R17  R21  R25
7  -4  5  31                      R14  R18  R22  R26

Here, the control number of this array is   11.02604   and if we take p = 10-7 again,

11.02604  ENTER^
E-7     XEQ "LS"   gives   -5.4  10-17  in X-register and in R00 = the determinant of the 4x4 matrix, which is of course 0   ( execution time = 21s )

and registers R11 thru R27 are now:

1  0   0.333333333     2.333333333
0  1  -0.666666667   -3.666666669
0  0      5.10-10              -2.10-9
0  0          0                    4.10-9

Thus, the system is equivalent to:

x + z /3  =  7/3              and there is an infinite set of solutions:
y - 2z /3 = -11/3           z may be chosen at random and x and y are determined by
0  =  0                 x = -z /3 + 7/3
0  =  0                 y = 2z /3 - 11/3

Note:    If the last number of the initial matrix is 32 instead of 31 the array is changed into:

1  0   0.333333333   0        i-e         x + z /3 = 0
0  1  -0.666666667   0                     y - 2z /3 = 0
0  0        5.10-10        0                                 0 = 0
0  0           0              1                                 0 = 1        and there is no solution at all !

Example3:   Invert Pascal's matrix:

1  1  1   1   1
1  2  3   4   5
1  3  6  10 15
1  4 10 20 35
1  5 15 35 70

This is equivalent to solve 5 linear systems simultaneously. If you begin at register R01, you have to store the 50 numbers

1  1  1   1  1    1  0  0  0  0                  R01  R06  R11  R16  R21  R26  R31  R36  R41  R46
1  2  3   4  5    0  1  0  0  0                  R02  R07  R12  R17  R22  R27  R32  R37  R42  R47
1  3  6  10 15  0  0  1  0  0         in      R03  R08  R13  R18  R23  R28  R33  R38  R43  R48        ( the right half is the unit matrix )
1  4 10 20 35  0  0  0  1  0                  R04  R09  R14  R19  R24  R29  R34  R39  R44  R49
1  5 15 35 70  0  0  0  0  1                  R05  R10  R15  R20  R25  R30  R35  R40  R45  R50

The control number of this rectangular array is 1.05005 and we can choose p = 0 because Pascal's matrix is non-singular:

1.05005  ENTER^
0      XEQ "LS"  yields  1 in X-register and in R00  ( the determinant of Pascal's matrices is always 1 )         ( execution time = 79s )

The array is now:

1  0  0  0  0   5  -10  10  -5   1
0  1  0  0  0 -10  30 -35  19 -4
0  0  1  0  0  10 -35  46 -27  6           ( the unit matrix is now on the left and the inverse matrix is on the right,
0  0  0  1  0  -5   19 -27  17 -4              in registers R21 thru R50 )
0  0  0  0  1   1   -4    6   -4   1

-This program can be used to invert a n x n matrix, but it requires 2n2 registers and thus, it cannot invert a 13x13 matrix.
-See "Matrix Operations for the HP-41" if you want to find the inverse of up to a 17x17 matrix.

b)  Program#2

-This variant below performs the same calculations as "LS"  but the coefficients are to be stored from register R01
-If flag F00 is clear the pivots smaller than E-7 ( line 60 ) is regarded as zero
-If flag F00 is set, only zero is regarded as zero
-Like "LS" , CF 01 = partial pivoting ,  SF 01 = no pivoting

-Line 60,  E-7  is a small number that identifies tiny elements.
-Line 108 is a three-byte GTO 01

 01  LBL "LS2"    02  X<>Y   03  .1   04  %   05  +   06  STO N   07  FRC   08  STO O   09  *   10  1   11  STO 00          12  ST+ O   13  LASTX   14  %   15  +   16  +   17  LBL 01   18  STO M   19  FS? 01   20  GTO 04   21  INT   22  RCL O   23  FRC 24  +   25  ENTER^   26  ENTER^   27  CLX   28  LBL 02   29  RCL IND Z   30  ABS   31  X>Y?   32  STO Z   33  X>Y?   34  +   35  RDN   36  ISG Z   37  GTO 02   38  RCL M   39  ENTER^   40  FRC   41  R^   42  INT   43  +   44  X=Y?   45  GTO 04   46  LBL 03 47  RCL IND X   48  X<> IND Z   49  STO IND Y   50  ISG Y   51  RDN   52  ISG Y   53  GTO 03   54  RCL 00   55  CHS   56  STO 00   57  LBL 04   58  CLX   59  FC? 00   60   E-7    61  RCL IND M   62  ST* 00   63  ABS   64  X<=Y?   65  GTO 09   66  RCL M   67  LASTX   68  LBL 05   69  ST/ IND Y 70  ISG Y   71  GTO 05   72  RCL O   73  STO P    74  LBL 06   75  RCL M   76  ENTER^   77  FRC   78  RCL P   79  INT   80  +   81  X=Y?   82  GTO 08   83  RCL IND X   84  SIGN   85  RDN   86  LBL 07   87  RCL IND Y   88  LASTX   89  *   90  ST- IND Y   91  ISG Y   92  RDN 93  ISG Y   94  GTO 07   95  LBL 08   96  ISG P   97  GTO 06   98  LBL 09          99  RCL N 100  ST+ O 101  X^2 102  1 103  RCL M  104  + 105  X>Y? 106  CLX 107  ISG X 108  GTO 01  109  RCL 00 110  CLA 111  END

( 168 bytes / SIZE n.m+1 )

 STACK INPUTS OUTPUTS Z / ( n.nnn )2 Y n / X m det A

n is the number of rows
m is the number of columns

Z-output is useful to retrieve n
If n <= m  L-output = n2.nm,nn

Example:   Of course, all the examples given after the "LS" listing can be solved ( but the first coefficient must be stored into R01 )  for instance:

2x + 3y + 5z + 4t = 39
-4x + 2y + z + 3t  = 15
3x  -  y + 2z + 3t = 19
5x + 7y - 3z + 2t = 18

-We store these 20 numbers:

2  3  5  4  39                  R01  R05  R09  R13  R17
-4  2  1  3  15        =        R02  R06  R10  R14  R18      respectively
3 -1  2  3  19                  R03  R07  R11  R15  R19
5  7 -3  2  18                  R04  R08  R12  R16  R20

-There are 4 rows and 5 columns,

CF 00   CF 01
4  ENTER^
5  XEQ "LS2"  >>>>  Det A = 840 = R00  ( in 30 seconds )

-Registers R01 thru R16 now contains the unit matrix and registers R17 thru R20 contains the solution  x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:

1  0  0  0  1
0  1  0  0  2              the solution is the last column
0  0  1  0  3              of the new matrix.
0  0  0  1  4

-"LS2" is slightly faster than "LS"
-When the program stops, R00 = det A

Note:   If you prefer to choose the "small" number that identifies the tiny elements:

-Replace register M by register Q
-Replace line 102 by ENTER^   SIGN
-Replace lines 58 to 60 by  RCL M
-Replace line 02 by STO M   X<> Z

-In this case, the inputs must be:

 STACK INPUTS OUTPUTS Z n ( n.nnn )2 Y m / X p det A

n = the number of rows
m = the number of columns
p = the "small" number you have chosen.

c)  Program#3

-In this third variant, the matrix is the right-part of the array so that, when the program stops,
the solution ( x1 , x2 , ............. , xn ) is in registers  R01  R02  ..............  Rnn
-Here, the attempt to diagonalization starts by the lower right corner of the matrix.
-Flags F00 & F01 play the same role as above.

-Line 113 is a three-byte GTO 01

 01  LBL "LS3"   02  X<>Y   03  .1   04  %   05  +   06  STO N   07  ST* Y   08  FRC   09  -   10  STO O   11  1   12  STO 00          13  LASTX   14  %   15  RCL Z   16  INT   17  +   18  LBL 01   19  STO M   20  FS? 01   21  GTO 04   22  INT   23  RCL O   24  FRC 25  +   26  ENTER^   27  ENTER^   28  CLX   29  LBL 02   30  RCL IND Z   31  ABS   32  X>Y?   33  STO Z   34  X>Y?   35  +   36  RDN   37  DSE Z   38  GTO 02   39  RCL M   40  ENTER^   41  FRC   42  R^   43  INT   44  +   45  X=Y?   46  GTO 04   47  LBL 03   48  RCL IND X 49  X<> IND Z   50  STO IND Y   51  DSE Y   52  RDN   53  DSE Y   54  GTO 03   55  RCL 00   56  CHS   57  STO 00   58  LBL 04   59  CLX   60  FC? 00   61   E-7   62  RCL IND M   63  ST* 00   64  ABS   65  X<=Y?   66  CLX   67  X=0?   68  STO 00   69  X=0?   70  GTO 09   71  RCL M   72  LASTX 73  LBL 05   74  ST/ IND Y   75  DSE Y   76  GTO 05   77  RCL O   78  STO P   79  LBL 06   80  RCL M   81  ENTER^   82  FRC   83  RCL P   84  INT   85  +   86  X=Y?   87  GTO 08   88  RCL IND X   89  SIGN   90  RDN   91  LBL 07   92  RCL IND Y   93  LASTX   94  *   95  ST- IND Y   96  DSE Y 97  RDN   98  DSE Y   99  GTO 07 100  LBL 08 101  DSE P 102  GTO 06        103  LBL 09 104  RCL N 105  ST- O 106  RCL O 107  RCL M 108  1 109  - 110  X<=Y? 111  CLX 112  DSE X 113  GTO 01 114  RCL 00 115  CLA 116  END

( 175 bytes / SIZE  n.m+1 )

 STACK INPUTS OUTPUTS T / n.nnn Z / / Y n / X m det A

n  =  number of rows
m = number of columns

T-output is useful to retrieve n

Example:

2x + 3y + 5z + 4t = 39                                       39 = 2x + 3y + 5z + 4t
-4x + 2y + z + 3t  = 15           is re-written          15 = -4x + 2y + z + 3t
3x  -  y + 2z + 3t = 19                                       19 = 3x  -  y + 2z + 3t
5x + 7y - 3z + 2t = 18                                       18 = 5x + 7y - 3z + 2t

and we store these 20 numbers:

39    2  3  5  4                        R01  R05  R09  R13  R17
15  -4  2  1  3         in            R02  R06  R10  R14  R18      respectively
19   3 -1  2  3                        R03  R07  R11  R15  R19
18   5  7 -3  2                        R04  R08  R12  R16  R20

-There are 4 rows and 5 columns,

CF 00   CF 01
4  ENTER^
5  XEQ "LS3"  >>>>  Det A = 840 = R00  ( in 30 seconds )

-Registers R05 thru R16 now contain the unit matrix and registers R01 thru R04 contain the solution  x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:

1  1  0  0  0
2  0  1  0  0                the solution is the first column
3  0  0  1  0                of the new matrix.
4  0  0  0  1

-When the program stops, R00 = det A
-If you have to invert a matrix, the inverse will be in registers R01 thru Rn2

-Moreover, the programs "LS-" & "LS+" below become shorter if they call "LS3" instead of "LS2"

d)  Underdetermined Systems

-"LS" & "LS2" can find the solution(s) of an underdetermined system A.X = B ( if any )
-But the following programs compute the Moore-Penrose solution given by the formula:   X = tA ( A.tA ) -1 B

d-1)  Program#1

-Store the coefficients of the system in column order , starting with register R01
-Flags F00 & F01 are used like with "LS2"
- "LS2" cf paragraph b) above
"TRN" transpose of a matrix                                     ( cf "Matrix Operations for the HP-41" )
"M*"  multiplication of 2 matrices
"TM*" multiplication of the transpose of a matrix by another one
"MCO" copying a matrix
are called as subroutines

-Lines 34 thru 44 may be replaced by

RCL 00        /              1           +
LASTX       +              +          REGMOVE
+                 RCL 00    E3
E6               X^2          /

-If you prefer to use "LS" instead of "LS2" replace lines 45 to 51 by

1               E2         /         XEQ "LS"
RCL 00      /           1         LASTX
ST+ Y        +          +         FRC
ST* Y       E3      E-7       ISG X

 01  LBL "LS-"  02  STO 00        03  STO Z 04  DSE X 05  X<>Y 06  ST* Y 07  ST* Z 08   E2 09  / 10  + 11   E3 12  / 13  1 14  + 15  X<>Y 16  2 17  + 18  XEQ "TRN" 19  ENTER^ 20  STO Z 21  1 22  - 23  RDN 24  STO IND T 25  LASTX 26  XEQ "TM*" 27  X<>Y 28  ENTER^ 29  X<> 00        30  DSE X 31  * 32  1 33  + 34  STO Y 35  RCL 00 36  +  37   E3  38  /  39  +  40  RCL 00 41  X^2 42  1 43  + 44  XEQ "MCO" 45  RCL 00 46  RCL 00 47  1  48  +  49  XEQ "LS2" 50  X<> Z  51  SQRT 52  INT 53  STO 01 54  X^2 55  STO Y 56  LASTX 57  + 58  1 59  + 60  RCL IND X 61  STO T 62  SIGN 63  - 64  E3 65  / 66  + 67  ISG X 68  RCL 01        69   E5 70  / 71  + 72  1 73  XEQ "M*"  74  END

( 120 bytes / SIZE 2n.m-n+2 )

 STACK INPUTS OUTPUTS Y n / X m 1.rrr,rr

n = number of rows
m = number of columns      ( n < m-1 )
1.rrr,rr  is the control number of the solution ( in fact r = m-1 )

Example:    Let's find the Moore-Penrose solution of the following system:

2x + 3y + 7z + 4t = 1                2  3   7   4  1           R01  R04  R07  R10  R13
3x + 2y - 5z + 8t = 4     store    3  2  -5  8  4   into   R02  R05  R08  R11  R14   respectively
4x + 5y + 6z + t  = 7                 4  5   6   1  7           R03  R06  R09  R12  R15

-There are 3 rows and 5 columns so:

3  ENTER^
5  XEQ "LS-"  >>>>   1.00404  ( in 49 seconds )

-The solution is in registers   R01 R02 R03 R04  namely:    x = 1.095088161 , y = 1.075776658 , z = -0.389378673 , t = -0.422963897

-When the program stops, R00 = det ( A.tA )   [ in this example, R00 = 128628 ]
-If R00 = 0 or is very "small" ( suggesting A.tA is singular ) , the results are probably meaningless.

d-2)  Program#2

-This variant also calls   "TRN"  "M*"  "TM*"  "MCO"   as subroutines    ( cf "Matrix Operations for the HP-41" )
but it uses "LS3" instead of "LS" or "LS2"

 01  LBL "LS3-"  02  STO 00 03   E-2 04  + 05   E3 06   / 07  RCL Y 08  ST* Y 09  + 10  RCL 00        11  RCL Z 12  * 13  1 14  ST+ Z 15  + 16  XEQ "TRN" 17  ENTER^ 18  STO 00 19  X<> Z 20  FRC 21  ISG X 22  INT 23  1 24  + 25  XEQ "TM*"  26  LASTX 27  ENTER^ 28  ENTER^ 29  DSE X 30  STO T 31  * 32  1 33  + 34  X<> 00 35  STO IND 00 36  RDN 37  XEQ "LS3" 38  R^ 39  INT 40  ENTER^ 41  ENTER^ 42  X^2 43  + 44  1 45  + 46  RCL IND X  47  ENTER^ 48  INT 49  DSE X 50  R^ 51  1 52  % 53  + 54   E3 55  / 56  1 57  + 58  X<>Y 59  XEQ "M*" 60  1 61  XEQ "MCO" 62  END

( 106 bytes / SIZE 2n.m-n+1 )

 STACK INPUTS OUTPUTS Y n / X m 1.rrr,rr

n = number of rows
m = number of columns      ( n < m-1 )
1.rrr,rr  is the control number of the solution ( in fact r = m-1 )

Example:    Let's find the Moore-Penrose solution of the same system:

1 = 2x + 3y + 7z + 4t                 1  2  3   7   4           R01  R04  R07  R10  R13
4 = 3x + 2y - 5z + 8t      store    4  3  2  -5  8   into   R02  R05  R08  R11  R14   respectively
7 = 4x + 5y + 6z + t                   7  4  5   6   1           R03  R06  R09  R12  R15

-There are 3 rows and 5 columns so:

3  ENTER^
5  XEQ "LS3-"  >>>>   1.00404  ( in 47 seconds )

-The solution is in registers   R01 R02 R03 R04  namely:    x = 1.095088162 , y = 1.075776659 , z = -0.389378674 , t = -0.422963896

-When the program stops, R00 = det ( A.tA )   [ in this example, R00 = 128628 ]
-If R00 = 0 ( suggesting A.tA is singular ) , the results are probably meaningless.

e)  Overdetermined Systems

-An overdetermined system  A.X = B  has often no solution at all!
-But the following programs calculate the least-squares solution:
-It minimizes || A.X - B ||

Formula:   X = ( tA.A ) -1 ( tA.B )

e-1)  Program#1

-Store the coefficients of the system in column order , starting with register R01
-Flags F00 & F01 are used like with "LS2"
- "LS2"  cf paragraph b) above
"TM*" multiplication of the transpose of a matrix by another one             ( cf "Matrix Operations for the HP-41" )
"MCO" copying a matrix
are called as subroutines

-If you prefer to use "LS" instead of "LS2" replace lines 25 to 33 by

RCL 00       E2         /         XEQ "LS"       INT
ENTER^      /           1         LASTX           X^2
DSE X        +          +         FRC                STO Y
ST* Y         E3       E-7      ISG X             LASTX

 01  LBL "LS+"  02  STO 00        03  RCL Y 04  ST* Y 05   E2 06  / 07  + 08  RCL X 09  RCL Z 10  - 11   E3 12  ST/ Z 13  / 14  1 15  ST+ Z 16  + 17  X<> Z 18  RCL 00 19  * 20  1 21  + 22  XEQ "TM*" 23  1 24  XEQ "MCO" 25  RCL 00  26  DSE X 27  RCL 00 28  XEQ "LS2"    29  X<> Z  30  INT  31  ENTER^ 32  ENTER^ 33  SQRT 34  + 35   E3 36  / 37  + 38  1 39  ST+ Y 40  XEQ "MCO" 41  END

( 78 bytes / SIZE n.m+m2 -m+1 )

 STACK INPUTS OUTPUTS Y n / X m 1.rrr,01

n = number of rows
m = number of columns      ( n >= m )
1.rrr,01 is the control number of the solution ( in fact r = m-1 )

Example:       The system:

5x + y + z  =   8             has no solution.
4x - y + 2z =  13
x + 2y - z =  -5
7x - 4y + 5z =  32
2x + 5y - 9z = -20

-We have to store these 20 coefficients into registers R01 thru R20 ( in column order )

5   1   1   8                       R01  R06  R11  R16
4  -1  2  13                      R02  R07  R12  R17          respectively
1   2  -1 -5          =          R03  R08  R13  R18
7  -4  5  32                      R04  R09  R14  R19
2   5 -9 -20                      R05  R10  R15  R20

-There are 5 rows and 4 columns so

5  ENTER^
4  XEQ "LS+"  >>>>   1.00301  ( in 54 seconds )

-The "least-squares solution" is in registers    R01 R02 R03  namely:   x = 2.071207430  ,  y = -3.201238385  ,  z = 0.904024772

-When the program stops, R00 = det ( tA.A )   [ In this example  R00 = 55233 ]
-If det ( tA.A )  = 0 ( suggesting tA.A is singular ) , the results are probably meaningless.

e-2)  Program#2

-This variant also calls  "TM*" and "MCO"   as subroutines    ( cf "Matrix Operations for the HP-41" )
but it uses "LS3" instead of "LS" or "LS2"

 01  LBL "LS3+"  02  STO 00        03   E3 04  / 05   E-5 06  + 07  1 08  + 09  RCL Y  10  * 11  ENTER^      12  FRC 13  1 14  ST+ Z 15  + 16  RCL 00 17  R^ 18  * 19  1 20  + 21  XEQ "TM*" 22  1 23  XEQ "MCO" 24  RCL 00 25  DSE X 26  RCL 00 27  XEQ "LS3" 28  R^ 29  FRC             30  ISG X 31  END

( 59 bytes / SIZE n.m+m2 -m+1 )

 STACK INPUTS OUTPUTS Y n det X m 1.rrr

n = number of rows
m = number of columns      ( n >= m )
1.rrr  is the control number of the solution ( in fact r = m-1 )

Example:      The system:

8  =  5x + y + z              still has no solution.
13 =  4x - y + 2z
-5 =   x + 2y - z
32 = 7x - 4y + 5z
-20 =  2x + 5y - 9z

-We store these 20 coefficients into registers R01 thru R20   ( in column order )

8    5   1   1                     R01  R06  R11  R16
13   4  -1  2                     R02  R07  R12  R17
-5   1   2  -1         =         R03  R08  R13  R18          respectively
32   7  -4  5                     R04  R09  R14  R19
-20   2   5  -9                    R05  R10  R15  R20

-There are 5 rows and 4 columns so

5  ENTER^
4  XEQ "LS3+"  >>>>   1.003  ( in 49 seconds )

-The "least-squares solution" is in registers    R01 R02 R03  namely:   x = 2.071207430  ,  y = -3.201238398  ,  z = 0.904024763

-When the program stops, Y = R00 = det ( tA.A )   [ In this example  R00 = 55233 ]
-If det ( tA.A )  = 0 or is very "small" ( suggesting tA.A is singular ) , the results are probably meaningless.

f)  Tridiagonal Systems

-The following program solves a  nxn  linear system  A.X = B  where  A = [ ai,j ]  is a tridiagonal matrix  i-e  ai,j = 0  if  | i - j | > 1
-In other words, only the main diagonal and its 2 nearest neighbors have non-zero elements.

-Gaussian elimination is used without pivoting, so the method can fail, for instance if  a11 = 0
but practically, most of the problems leading to tridiagonal systems do not have these pitfalls.

-The components of the 3 diagonals are to be stored in column order into contiguous registers, starting with a register Rbb of your own choosing ( with bb  > 00 )
followed by the  bi 's ( the n elements of B )

-The determinant  det A  is also computed and stored in R00.

Data Registers:    R00 = det A  at the end                                              ( Registers Rbb thru Ree are to be initialized before executing "3DLS" )

•  Rbb     = a11    •  Rbb+2 = a12                                                                                                                          •  R3n+bb-2 = Ree-n+1 = b1
•  Rbb+1 = a21   •  Rbb+3 = a22    •  Rbb+5 = a23                                                                                                •  R3n+bb-1 = Ree-n+2 = b2
•  Rbb+4 = a32    •  Rbb+6 = a33         ............
•  Rbb+7 = a43         ............                                                                                      .....................
............

•  R3n+bb-10 = an-3,n-2
•  R3n+bb-9 = an-2,n-2      •  R3n+bb-7 = an-2,n-1
•  R3n+bb-8 = an-1,n-2      •  R3n+bb-6 = an-1,n-1    •  R3n+bb-4 = an-1,n
•  R3n+bb-5 = an,n-1      •  R3n+bb-3 = an,n            •  R4n+bb-3 = Ree = bn

>>>>  When the program stops, the solutions  x1 , .... , xn  have replaced  b1 , .... , bn  in registers  R3n+bb-2  thru  R4n+bb-3

Flags: /
Subroutines: /

 01  LBL "3DLS" 02  INT 03  LASTX 04  FRC 05   E3 06  * 07  RCL X 08  3 09  + 10  R^ 11  - 12  4 13  / 14  STO M 15  - 16  1 17  STO 00 18  + 19  STO N 20   E3 21  / 22  + 23  STO O 24  LBL 01 25  RCL O 26  RCL IND X 27  ST* 00 28  ST/ IND N 29  ISG O 30  ISG O 31  FS? 30 32  GTO 02 33  ST/ IND O 34  CLX 35  SIGN 36  + 37  RCL IND N 38  RCL IND Y 39  * 40  ISG N 41  CLX 42  ST- IND N 43  RCL IND O 44  LASTX 45  * 46  ISG O 47  ST- IND O 48  GTO 01 49  LBL 02 50  CLX 51  SIGN 52  - 53  INT 54  3 E-5 55  + 56  RCL M 57  DSE X 58  STO O 59  X<>Y 60  LBL 03 61  RCL IND X 62  RCL IND N 63  * 64  DSE N 65  ST- IND N  66  DSE Y 67  X<>Y 68  DSE O 69  GTO 03 70  RCL N 71  RCL M 72  + 73  DSE X 74   E3 75  ST/ Y 76  RDN 77  ST+ N 78  X<> N 79  CLA 80  END

( 132bytes / SIZE 4n+d-2 )

 STACK INPUTS OUTPUTS X (bbb.eee)system (bbb.eee)solution L / n

(bbb.eee)system  is the control number of the system  ,  eee =  4n + bbb - 3  , bb > 00 ,  n = number of unknowns ,  n > 1
(bbb.eee)solution is the control number of the solution

Example:                                                                                   With  bb = 01,  store                                            into

2.x + 5.y                                     = 2                                       2   5                           2                 R01   R03                                           R17
3.x + 7.y + 4.z                            = 4                                       3   7   4                      4                 R02   R04   R06                                  R18
y + 3.z + 7.t                    = 7                                            1   3   7                 7                          R05   R07   R09                         R19
2.z + 4.t + 6.u           = 1                                                 2   4   6            1                                    R08   R10   R12                R20
8.t +  u   + 7.v  = 5                                                      8   1   7       5                                             R11   R13   R15       R21
9.u + 4.v  = 6                                                           9   4       6                                                       R14   R16       R22

-Then,   1.022  XEQ "3DLS"  >>>>   17.022    ( in 9 seconds )

R17 = x = -16.47810408            R20 = t = -0.480422463
R18 = y =     6.991241632          R21 = u =  0.112313241
R19 = z =     1.123905204          R22 = v =  1.247295209

-We have also  R00 = det A = 3882

Notes:

-Synthetic registers M N O may be replaced by R01 R02 R03 but in this case, choose bb > 03
-The execution time is approximately proportional to n.
-With bb = 01, this program can solve a tridiagonal system of 80 equations in 80 unknowns in about  118 seconds,  provided it is executed from extended memory.
-There is a risk of OUT OF RANGE  line 27 if the determinant exceeds 9.999999999 E99 - especially for large n-values. Set flag F24 in this case.
-If you don't want to calculate det A , delete lines 27 and 17 , register R00 will be unused.

-"5DLS" hereunder solves a  nxn  linear system  A.X = B  where  A = [ ai,j ]  is a pentadiagonal matrix  i-e  ai,j = 0  if  | i - j | > 2
-Only the main diagonal and its 4 nearest neighbors have non-zero elements.

-Gaussian elimination is used without pivoting.

-The components of the 5 diagonals are to be stored in column order into contiguous registers, starting with a register Rbb of your own choosing ( with bb  > 00 )
followed by the  bi 's ( the n elements of B )

-The determinant  det A  is computed and stored in R00.

Data Registers:    R00 = det A  at the end                                                 ( Registers Rbb thru Ree are to be initialized before executing "5DLS" )

•  Rbb     = a11    •  Rbb+3 = a12   •  Rbb+7 = a13                                                                                                     •  R5n+bb-6 = Ree-n+1 = b1
•  Rbb+1 = a21   •  Rbb+4 = a22    •  Rbb+8 = a23     •  Rbb+12 = a24                                                                       •  R5n+bb-5 = Ree-n+2 = b2
•  Rbb+2 = a32   •  Rbb+5 = a32    •  Rbb+9 = a33     •  Rbb+13 = a34
•  Rbb+6 = a42    •  Rbb+10 = a43   •  Rbb+14 = a44
•  Rbb+11 = a53   •  Rbb+15 = a54
•  Rbb+16 = a64

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

•  R5n+bb-23 = an-5,n-3
•  R5n+bb-22 = an-4,n-3     •  R5n+bb-18 = an-4,n-2
•  R5n+bb-21 = an-3,n-3     •  R5n+bb-17 = an-3,n-2   •  R5n+bb-13 = an-3,n-1
•  R5n+bb-20 = an-2,n-3     •  R5n+bb-16 = an-2,n-2  •  R5n+bb-12 = an-2,n-1    •  R5n+bb-9 = an-2,n
•  R5n+bb-19 = an-1,n-3     •  R5n+bb-15 = an-1,n-2  •  R5n+bb-11 = an-1,n-1    •  R5n+bb-8 = an-1,n
•  R5n+bb-14 = an,n-2     •  R5n+bb-10 = an,n-1      •  R5n+bb-7 = an,n              •  R6n+bb-7 = Ree = bn

>>>>  When the program stops, the solutions  x1 , .... , xn  have replaced  b1 , .... , bn  in registers  R5n+bb-6  thru  R6n+bb-7

Flags:  F06 & F07  are cleared by this routine.
Subroutines: /

-Line 96 is a three-byte  GTO 01

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

( 265 bytes / SIZE 6n+bb-6 )

 STACK INPUTS OUTPUTS X (bbb.eee)system (bbb.eee)solution L / n

(bbb.eee)system  is the control number of the system  ,  eee =  6n + bbb - 7  ,  bb > 00 ,  n = number of unknowns ,  n > 2
(bbb.eee)solution is the control number of the solution

Example:

7 x + 3 y + 4 z                                          =  1
x + 8 y + 6 z +  t                                    =  2
3 x + 2 y + 9 z + 2 t + 3 u                         =  3
4 y + 4 z + 8 t + 5 u + 2 v                =  4
2 z + 3 t + 9 u + 3 v +  w        =  5
2 t + 3 u + 7 v + 2 w      =  6
u + 6 v + 8 w      =  7

-With  bb = 07  store these 36 numbers

7   3   4                           1                             R07   R10   R14                                             R36
1   8   6   1                      2                             R08   R11   R15   R19                                    R37
3   2   9   2   3                 3                             R09   R12   R16   R20   R24                           R38
4   4   8   5   2            4               in                     R13   R17   R21   R25   R29                 R39                respectively
2   3   9   3   1       5                                                R18   R22   R26   R30   R33        R40
2   3   7   2       6                                                          R23   R27   R31   R34        R41
1   6   8       7                                                                   R28   R32   R35        R42

7.042  XEQ "5DLS"  >>>>     36.042    ( in 23 seconds )

R36 = x = -0.023088805         R39 = t = 0.038788731       R42 = w = 0.364890723
R37 = y =   0.069105035         R40 = u = 0.235429731
R38 = z =   0.238576632         R41 = v = 0.640907414

-We have also  R00 = det A = 607264

Notes:

-Synthetic registers M N O P Q may be replaced by R01 R02 R03 R04 R05 but in this case, choose bb > 05
-"5DLS" cannot solve a 2x2 linear system: there must be at least 5 "diagonals".
-The execution time is roughly proportional to n.
-This program can solve a pentadiagonal system of 54 equations in 54 unknowns ( in 197 seconds ),  provided it is executed from extended memory.
-There is a risk of OUT OF RANGE  line 37 or 102 if the determinant exceeds 9.999999999 E99 - especially for large n-values. Set flag F24 in this case.
-If you don't want to calculate det A , delete lines 102 , 37 and 24 , register R00 will be unused.

2°) Non-Linear Systems

a) One Equation in one Unknown:   f(x) = 0

a-1) Secant Method

-Although it's obviously not a system of equations, it is logical to begin with this short routine.

Data Registers:           •  R00 = function name                                ( Register R00 is to be initialized before executing "SOLVE" )

R01 = x , R02 = x' , R03 = f(x)
Flags: /
Subroutine:  A program which computes f(x) assuming x is in X-register upon entry.

 01  LBL "SOLVE" 02  STO 01 03  X<>Y 04  STO 02 05  XEQ IND 00 06  STO 03 07  LBL 01 08  VIEW 01 09  RCL 01 10  XEQ IND 00 11  ENTER^ 12  ENTER^ 13  X<> 03 14  - 15  X#0? 16  / 17  RCL 02 18  RCL 01 19  STO 02 20  STO T 21  - 22  * 23  + 24  STO 01 25  X#Y? 26  GTO 01 27  END

( 43 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS T / x Z / x Y x'0 x X x0 x

where  x0  and  x'0  are 2 initial guesses with  x0 # x'0

Example:   Find a solution near  x = 2 of  x3 - 4x + 1 = 0

LBL"T"
ENTER^
X^2
4
-
*
1
+
RTN      "T"  ASTO 00  and if we choose  2 and 3 as initial guesses:

2  ENTER^
3  XEQ "SOLVE"  >>>>  the successive x-values are displayed and eventually:

x = 1.860805853 = X-register = R01 = R02  ( f(x) ~ 10 -10 = R03 )

-The other solutions are  x' = 0.254101688  and  x" = -2.114907542   ( the last decimal should be a 1 )

Note:

-Always check f(x) in R03 because the iterations will also stop if  f(x) = f(x') , thus producing a possibly wrong result!

-The next step in interpolation methods is to use a parabola instead of a straight-line.
-Three starting-values are needed:   x1 , x2 , x3   The next approximation x4 is computed by the formula:

x4 = x3 + 2.C / [ -B ± ( B2 - 4 A.C )1/2 ]     ( where the sign is chosen to produce the largest denominator )   with:

A = ( x2 - x1 ) f(x3) + ( x1 - x3 ) f(x2) + ( x3 - x2 ) f(x1)
B = ( x2 - x1 ) ( 2 x3 - x2 - x1 ) f(x3) - ( x1 - x3 )2 f(x2) + ( x3 - x2 )2 f(x1)
C = ( x3 - x2 ) ( x2 - x1 ) ( x3 - x1 ) f(x1)

-If, however, the discriminant is negative, x4 is computed by  x4 = x3 - B/(2A)
-This facilitates the determination of double-roots.

-And sometimes , the algorithm leads to a local extremum.

Data Registers:           •  R00 = function name                       ( Register R00 is to be initialized before executing "SOLVE2" )

R01 = x3      R04 = f(x3)       R07 thru R10: temp
R02 = x2      R05 = f(x2)
R03 = x1      R06 = f(x1)
Flags: /
Subroutine:  A program which computes f(x) assuming x is in X-register upon entry.

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

( 113 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS Z x1 / Y x2 f(x) X x3 x

Example:    f(x) = 12 x3 - 44 x2 - 5 x + 100

LBL "T"       12       -       -        +
ENTER^      *        *       *      RTN
STO Z         44       5     100

"T"  ASTO 00

-If the starting-values are  -3  -2  -1

-3  ENTER^
-2  ENTER^
-1  XEQ "SOLVE2"  >>>>  the successive approximations are displayed and finally:

x = -1.333333333   RDN   f(x) = 0.000000050     ( in 17 seconds )

-If the starting-values are 1  2  3

1  ENTER^
2  ENTER^
3     R/S    >>>>  the successive approximations are displayed and:

x = 2.499981841   RDN    f(x) = 0     ( in 22 seconds )

-The double-root is actually  x = 2.5  and - as usual in such cases - the result is not very accurate.

Note:

-The rate of convergence is better than that of the secant method.
-Actually, it's the same as that of Newton's method: a quadratic convergence.

a-3) Ridders' Method

-If a function is defined for a < x < b , interpolation methods may produce an approximation outside the interval [ a , b ] thus causing a DATA ERROR.
-Ridders' method avoids this pitfall: if a root is bracketed beween a and b , the successive approximations remain bracketed.
-Starting with 2 initial guesses x1 & x2 so that  f(x1) . f(x2) < 0 ,  the formulas are:

x3 = ( x1 + x2)/2
x4 = x3 + ( x3 - x1 )  f(x3)  { [ f(x3) ]2 -  f(x1) f(x2) } -1/2  sign [  f(x1) -  f(x2) ]

Data Registers:          •  R00 = function name                                    ( Register R00 is to be initialized before executing "RIDDERS" )

R01 = x2          R03 = x1         R05 = x3         R07 = x4
R02 = f(x2)      R04 = f(x1)      R06 = f(x3)      R08 = E-8
Flags: /
Subroutine:  A program that computes f(x) assuming x is in X-register upon entry.

-Line 88  may be replaced by the M-Code routine  X#Y??  to avoid a possible infinite loop  ( cf "A few M-Code routines for the HP-41" )

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

( 127 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS Y x1 f(x) ~ 0 X x2 x

Example:    Find a solution of the equation  ( 1 - x2 ) -1/2 - 4 = 0  between  0  and  0.99

LBL "T"    LASTX    1/X       RTN
X^2          -               4
SIGN       SQRT       -

'T'  ASTO 00

0     ENTER^
0.99  XEQ "RIDDERS"  >>>>  the successive x-approximations are displayed and finally:

x = 0.968245837   RDN   f(x) ~ 0.000000003         --- execution time = 18 seconds ---

Note:

-The starting-values should verify  f(x1) f(x2) < 0
-Otherwise, you'll probably get a DATA ERROR message line 38 or 39
-However, it may also work ... sometimes!

a-4)  Discrete Data ( 5 points )

-When the function is only known by tabulated values, one can use "SOLVE" & "LAGR" ( cf "Lagrange Interpolation for the HP-41" )  to solve  f(x) = 0
-The following routine is much less general:
-We assume that the abscissas are equally spaced  xi+1 - xi = h  and that there are 5 data points   y(x1) = y1 , ........... , y(x5) = y5

-In this case, the Lagrange interpolation polynomial may be written  y = y3 + a n + b n2 + c n3 + d n4   with  n = ( x - x3 )/h
and  f(x) = 0  is re-written   n = - (  y3 + b n2 + c n3 + d n4 ) / a
-So, this method may fail! For instance when a = 0

Data Registers:              R00 = h                             ( Registers R00 thru R05 are to be initialized before executing "SLVD" )

•  R01 = y1
•  R02 = y2
•  R03 = y3                            R06 to R09: temp      R10 = x3
•  R04 = y4
•  R05 = y5
Flags: /
Subroutines: /

 01  LBL "SLVD" 02  STO 10        03  X<>Y 04  STO 00 05  RCL 02 06  RCL 04 07  - 08  STO 06 09  LASTX 10  RCL 02 11  + 12  RCL 03 13  ST+ X 14  - 15  12 16  ST* 06 17  * 18  RCL 01        19  RCL 05 20  + 21  RCL Y 22  3 23  / 24  - 25  RCL 03 26  ST+ X 27  - 28  STO 09 29  - 30  STO 07 31  RCL 05 32  RCL 06        33  6 34  / 35  + 36  RCL 01 37  - 38  ST+ X 39  ST+ 06 40  STO 08 41  CLX 42  LBL 01 43  ENTER^ 44  ENTER^ 45  VIEW X 46  RCL 09        47  * 48  RCL 08 49  + 50  * 51  RCL 07 52  + 53  * 54  * 55  RCL 03 56  24 57  * 58  + 59  RCL 06        60  / 61  X#Y? 62  GTO 01 63  RCL 00 64  * 65  RCL 10 66  + 67  CLD 68  END

( 87 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS Y h / X x3 x

Example:    Solve f(x) = 0  from the following data

 x -3 -1 1 3 5 y -7 -3 -1 4 10

-7   STO 01
-3   STO 02
-1   STO 03
4    STO 04
10   STO 05

2   ENTER^
1   XEQ "SLVD"  the successive n-values are displayed ( the first one is always 0 cf line 41 ) and eventually, it yields:

x = 1.534588700   --- execution time = 17 seconds ---

-With the starting-values 1 and 2 , "SOLVE" & "LAGR" produce the same result in 1mn51s

Notes:

-To avoid an infinite loop, you can replace line 61 by the M-code routine  X#Y??  ( cf "A few M-Code Routines for the HP-41" )

-Or replace lines 61-62 by  X<>Y  RND  X<>Y  RND  X#Y?  GTO 01  LASTX
and replace lines 41-42 by  CLX  +  LBL 01  LASTX

-The accuracy will be controlled by the display setting  (  FIX 9 gives the same result in the example above )

b) 2 Equations in 2 Unknowns:   f(x,y) = 0 ; g(x,y) = 0

- "SXY" uses the Newton's iterative method - more exactly a generalization of the secant method to approximate the partial derivatives.
Registers R00 thru R11 are used.
- It requires 2 initial guesses ( x , y ) and ( x' , y' ) which are to be stored in R01 , R02 and R03 , R04 respectively.
You must choose x # x'  and y # y'   ( very important )

- You also have to load a subroutine that takes y in Y-register and x in X-register
and calculates f(x;y) in Y-register and g(x;y) in X-register.

>>>   Thus, the stack has to be changed from:                  Y = y      into       Y' = f(x,y)             by your subroutine.
X = x                  X' = g(x,y)
R11 and greater are available for this job.
- Then, you store the name of this subroutine into R00   ( global label of 6 characters maximum )
and XEQ "SXY"

- The successive x-values are displayed  ( line 03 ) and when the program stops,

x is in X-register and in R01
y is in Y-register and in R02
Z-register contains | f(x,y) | + | g(x,y) |

f(x,y) is in R05
g(x,y) is in R06

- To find another solution, re-initialize R01 thru R04 and R/S

Warning:

-The program stops when the approximate Jacobian determinant equals zero.
-This happens when   x = x' or y = y' but it may also happen by a stroke of bad luck,
for instance if x converges much more quickly than y to the solution.

>>>  That's why it's always wise to check the value of  | f | +| g |

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

( 95 bytes / SIZE 012 )

 STACK INPUTS OUTPUTS Z / | f | + | g | Y / y X / x

Example:   Find x and y such that   x.y = 7 and  x2 + y4 = 30

First, we rewrite the system:

x.y -  7  = 0
x2 + y4 - 30 = 0

The subroutine may be ( let's call it "T" ):

01  LBL "T"
02  RCL Y
03  RCL Y
04  *
05  7
06  -
07  X<>Y
08  X^2
09  R^
10  X^2
11  X^2
12  +
13  30
14  -
15  RTN

"T"  ASTO 00     and  with ( 2 ; 2 ) ( 3 ; 3 ) as initial guesses
2  STO 01  STO 02
3  STO 03  STO 04
XEQ "SXY"              the successive x-values are displayed:

2.000000000
3.458333333
3.369648334
3.368074217
3.368200504
3.368200265
3.368200265

and the program stops with

X = 3.368200265 = x                              ( execution time = 26s )
Y = 2.078261222 = y
Z = 0.000000011 = | f(x;y) | + | g(x;y) |    ( generally a small number if it's a "true" solution )

c) 3 Equations in 3 Unknowns:   f(x,y,z) = 0 ; g(x,y,z) = 0 ; h(x,y,z) = 0

- The same method leads to the following program for solving a system of 3 non-linear equations.
Registers R00 thru R20 are used.
- Two initial guesses  ( x , y , z )  and  ( x' , y' , z' ) are to be stored into  R01 , R02 , R03 and R04 , R05 , R06 respectively   ( x#x' , y#y' , z#z' )
- You write a subroutine that takes   x in X-register ; y in Y-register ; z in Z-register
and  computes  f(x,y,z) in Z-register ,  g(x,y,z) in  Y-register  ,  h(x,y,z)  in X-register

The stack has to be change from:

Z = z                Z' = f(x,y,z)
Y = y     to       Y'= g(x,y,z)           ( registers R16 and greater are available for that purpose )
X = x               X'= h(x,y,z)

- Then, you store the name of this subroutine into R00 and XEQ "SXYZ"
The successive x-values are displayed and when the program stops,

x = X-register = R01
y = Y-register = R02
z = Z-register = R03

and   T-register contains | f(x,y,z) | + | g(x,y,z) | + | h(x,y,z) |

f(x,y,z) = R07
g(x,y,z) = R08
h(x,y,z) = R09

>> The recommendations for "SXY" still apply for "SXYZ"

Note:   The REGSWAP function is used in the following listing, but if you don't have an X-Functions module, you can:

a) delete line 91
b) replace lines 77 to 79 by  RCL 13   X<>16  STO 13  RCL 14  X<>17  STO 14  RCL 15  X<>18  STO 15
c) replace lines  65 to 67 by  RCL 10   X<>16  STO 10  RCL 11  X<>17  STO 11  RCL 12  X<>18  STO 12
d) replace lines  52 to 53 by  RCL 07   X<>16  STO 07  RCL 08  X<>17  STO 08  RCL 09  X<>18  STO 09
e) delete line 49

In this case, SIZE 020 is enough.

-Line 89 is a three-byte  GTO 01

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

( 189 bytes / SIZE 021 )

 STACK INPUTS OUTPUTS T / | f | + | g | + | h | Z / z Y / y X / x

Example:   Find a solution of the system:

x.y2 - z/y = 0
x - y - z   = 0
ln x + y.z = 0

The subroutines for calculating  f , g , h  may be:

01  LBL "T"
02  STO 16
03  X<>Y
04  STO 17
05  X^2
06  *
07  X<>Y
08  STO 18
09  RCL 17
10  /
11  -
12  RCL 16
13  RCL 17
14  -
15  RCL 18
16  -
17  RCL 17
18  LASTX
19  *
20  RCL 16
21  LN
22  +
23  RTN

"T"  ASTO 00                                  and, if you choose  ( 2,2,2 ) and (1,1,1) as initial approximations:

2  STO 01  STO 02  STO 03
1  STO 04  STO 05  STO 06
XEQ "SXYZ"

The successive x-values are displayed:

2.000000000
1.742625585
0.896281748
0.827584649
0.865942155
0.865417190
0.865408830
0.865408832
0.865408832

Finally, it yields:

x = 0.865408832 in X-register
y = 0.639295476 in Y-register
z = 0.226113356 in Z-register
| f | + | g | + | h | = 10-10        in T-register       ( the whole execution time is about 1mn37s )

More precisely,

f(x,y,z) = -10-10  in R07
g(x,y,z) =   0       in R08
h(x,y,z) =   0       in R09

>> Note that if you choose:

1 STO 01 STO 02 STO 03
2 STO 04 STO 05 STO 06  as initial guesses,

the program stops after only 2 iterations and returns x = 1 ; y = 0.777777778 ; z = 0.222222222
but you realize it's a completely wrong result because  T-register contains  0.492... ( not really very small! )

>> Good guesses are not always easy to find but in any case, always check the values of   f , g , h !

d) N Equations in N Unknowns:   f1( x1,... ,xn )  = 0 , .......... ,  fn( x1, ... ,xn ) = 0

- Now, we deal with the general case of a system of nxn non-linear equations.
These programs call "LS" , "LS2" or "LS3" as a subroutine ( synthetic version required ).
- Once again, the same ( quasi- ) Newton's method is used:
Each iteration solves a linear system of n equations in n unknowns and that's why "LS" is needed.
- Unlike "SXY" and "SXYZ"  it's of course impossible to take the n variables from the stack and calculate the n functions in the stack if n > 4,
therefore you'll have to key in n different subroutines for computing the fi in the X-register with x1 in R01 , ... , xn in Rnn
- Synthetic registers M N O  and data registers R00 thru Rn2+4n are used.

- The successive x1-values are displayed ( line 04 ) and when the program stops, | f1 | + .... + | fn | is in X-register
and the solution  ( x1 , .... , xn )  is  in  ( R01 , ..... , Rnn )

d-1)  Program#1

Data registers:        •  R00 = n                             ( Registers R00 thru R3n are to be initialized before executing "NLS" )

•  R01 = x1   •   Rn+1 = x'1      •  R2n+1 = f1  name
•  R02 = x2   •   Rn+2 = x'2      •  R2n+2 = f2  name                    R3n+1 thru Rn2+4n  contain the coefficients
.....................................................................                          of the linear system that "LS" solves at every iteration.

•  Rnn = xn   •     R2n = x'n      •  R3n  =  fn  name

where   (  x1 , .... , xn ) and ( x'1 , .... , x'n ) are the 2 initial guesses which must verify      xi # x'i   for  i = 1 , 2 , ... , n

Flag:   F01
Subroutines:     "LS"  ( synthetic version )  and  n programs that take  x1 in R01 , ..... , xn  in Rnn
and compute  fi ( x1 , .... , xn )  in the X-register   i = 1 , ..... , n

-Replace line 02 by  SF 01  if you want to solve the successive linear systems without pivoting:
it's of course faster but also more risky! ( it saves about 41 seconds in the example below )
-Line 107 is a three-byte GTO 01

 01  LBL "NLS"   02  CF 01    03  LBL 01    04  VIEW 01   05  4.00301   06  RCL 00   07  ST+ Y   08  *   09  STO M   10  3.002   11  LASTX   12  *   13  STO N   14  LBL 02   15  RCL IND N   16  XEQ IND X   17  LBL 03   18  STO IND M   19  DSE M   20  GTO 03   21  RCL 00   22  ENTER^   23  X^2   24  +   25  DSE X 26  ST+ M   27  DSE N   28  GTO 02   29  3   30  RCL 00   31  ST+ Y   32  *   33  STO M   34  3.002   35  LASTX   36  *   37  STO N   38  LASTX   39  .1   40  %   41  +   42  -   43  STO O   44  LBL 04   45  RCL O   46  RCL 00   47  -   48  RCL IND X   49  X<> IND O   50  STO IND Y 51  LBL 05   52  RCL IND N   53  XEQ IND X   54  ST- IND M   55  DSE M   56  DSE N   57  GTO 05   58  RCL O   59  RCL 00   60  ST+ N   61  -   62  RCL IND X   63  X<> IND O   64  STO IND Y   65  DSE O   66  GTO 04   67  2.01   68  RCL 00   69  ST+ Y   70  *   71   E3   72  /   73  RCL N   74  +   75  1 76  +   77  0   78  XEQ "LS"    79  LASTX    80  FRC    81  ISG X    82  INT   83  STO 00   84  STO M   85  ST+ X   86  STO N   87  ST+ X   88  RCL 00   89  X^2   90  +   91  STO O   92  X<>Y   93  X=0?   94  GTO 07   95  LBL 06   96  RCL IND M   97  ENTER^   98  X<> IND N   99  - 100  RCL IND O 101  * 102  ST- IND M 103  DSE O 104  DSE N 105  DSE M 106  GTO 06 107  GTO 01  108  LBL 07 109  CLA 110  3.002 111  RCL 00 112  * 113  STO N 114  LBL 08 115  RCL IND N 116  XEQ IND X 117  ABS 118  ST+ M 119  DSE N 120  GTO 08 121  X<> M 122  CLA 123  CLD 124  END

( 219 bytes / SIZE n2+4n+1 )

 STACK INPUTS OUTPUTS X / Sum | fi |

Example:    Find a solution near ( 4 , 1 , 3 , 6 ) of the system:

x1 + x2 + x3 + x4 - 16 = 0
x1 x2 x3 - 3 x4 = 0
4x12 - x2 x3 x4 - 40 = 0
x1 x2 x3  x4 - 140 = 0

-Let's write the 4 subroutines, "F1" "F2" "F3" "F4"  ( any global labels, 6 characters maximum ):

01  LBL "F1"
02  RCL 01
03  RCL 02
04  RCL 03
05  RCL 04
06  +
07  +
08  +
09  16
10  -
11  RTN
12  LBL "F2"
13  RCL 01
14  RCL 02
15  RCL 03
16  *
17  *
18  RCL 04
19  3
20  *
21  -
22  RTN
23  LBL "F3"
24  RCL 01
25  ST+ X
26  X^2
27  RCL 02
28  RCL 03
29  RCL 04
30  *
31  *
32  -
33  40
34  -
35  RTN
36  LBL "F4"
37  RCL 01
38  RCL 02
39  RCL 03
40  RCL 04
41  *
42  *
43  *
44  140
45  -
46  RTN

SIZE 033     ( or greater )
4  STO 00    ( the number of unknowns )    and  if     ( 4 , 1 , 3 , 6 )  and  ( 4.1 , 1.1 , 3.1 , 6.1 )  are your 2 initial guesses:

4  STO 01       4.1   STO 05     F1  ASTO 09
1  STO 02       1.1   STO 06     F2  ASTO 10
3  STO 03       3.1   STO 07     F3  ASTO 11
6  STO 04       6.1   STO 08     F4  ASTO 12

XEQ "NLS"    the successive x1-approximations are displayed:

4.000000000
4.298102981      ( each iteration lasts about 58 seconds )
4.266193620
4.266545269
4.266540470
4.266540475
4.266540475

and the program stops with  | f1 | + | f2 | + | f3 | + | f4 | = 10-8  in X-register  and the solution in registers  R01 , R02 , R03 , R04

Thus,   x1 =  4.266540475 , x2 =  1.353632236 ,  x3 = 3.548526779 ,  x4 = 6.831300511      ( the whole execution time is about 6mn33s )

>>>>  This system has several solutions. Another one is  ( 4.266540475 , -0.2552286465 , 18.81998868 , -6.831300511 )

Notes:

-I've tested this program to solve a 7x7 non-linear system of ( relatively ) simple equations: every iteration requires about 3 minutes ( see the exercise below ).
-Although I think it's rare, it might happen that this program runs for ever because of successive x-values that are very close together:
-In order to avoid this infinite loop:

Add  E-8 ( or another "small" number )  X<Y?  after line 106
Add  ABS  +  after line 102

-This modification also reduces execution time:
In the example above, the last ( unuseful ) iteration is avoided.

-There is another way to avoid the unuseful iteration ( when  xi = x'i for some i )

Replace lines 97 to 107 by

ENTER^  ENTER^  ENTER^  X<> IND N  -  RCL IND O  *  -  STO IND M  X=Y?  SF 01  DSE O  DSE N  DSE M  GTO 06  FC?C 01  GTO 01

-The same warnings mentioned for "SXY" and "SXYZ" still apply to "NLS"

*** If you want to use "LS2" instead of "LS", it's quite possible, but a little more complicated. For instance:

Replace lines 67 to 91 by

CLX  RCL 00  X^2  E3  /  +  E3  /  RCL M  +  1  +  REGSWAP  RCL 00  RCL 00  LASTX  +  XEQ "LS2"
RCL Z  INT  ENTER^  SQRT  STO 00  STO M  STO N  ST+ N  +  STO O  E3  /  ISG X  RCL 00  3  *
ST+ O  +  E3  /  1  +  REGSWAP

and add   SF 00   after line 01

-Actually, this variant is slightly faster: in the example above, execution time = 6m22s instead of 6m33s

Exercise:     Find a solution of the following system of 7 equations in 7 unknowns:

x3 + y2z + t.u - v2 - w2 = 0
x2y - z.t.u2 + x.v - w3 = 0
x + y + z + t - u - v - w = 0
x3 - y.z.t + t.u.v - w2 = 0
x.y4 - 2y.z3 - t.u.v2w = 0
x + y.z + t.u - v.w2 = 0
x.y - y.z.t.u.v + w - 1 = 0

Answer:    With the initial guesses   x = y = z = t = u = v = w = 1  &  x' = y' = z' = t' = u' = v' = w' = 2  ,  "NLS" yields after about 27 minutes:

x = 1.200271274         u = 1.280018254
y = 1.548046396         v = 1.698155509
z = 1.011876841         w = 1.463999879
t = 0.681979132

-With these values,  Sum i=1,....,7  | fi |  =  4 10 -9        ( Simply press XEQ 07 if you ever lose this result )

d-2)  Program#2

-Use the following variant if you prefer "LS3" to "LS" or "LS2"

Data registers:        •  R00 = n                             ( Registers R00 thru R3n are to be initialized before executing "NLS3" )

•  R01 = x1   •   Rn+1 = x'1      •  R2n+1 = f1  name
•  R02 = x2   •   Rn+2 = x'2      •  R2n+2 = f2  name
.....................................................................                                       R3n+1 to Rn2+4n:  temp

•  Rnn = xn   •     R2n = x'n      •  R3n  =  fn  name

where   (  x1 , .... , xn ) and ( x'1 , .... , x'n ) are the 2 initial guesses which must verify      xi # x'i   for  i = 1 , 2 , ... , n

Flags:  F00 - F01
Subroutines:     "LS3"  and  n routines that take  x1 in R01 , ..... , xn  in Rnn  and compute  fi ( x1 , .... , xn )  in X-register   i = 1 , ..... , n

-Line 127 is a three-byte  GTO 01

 01  LBL "NLS3"   02  CF 01   03  LBL 01   04  VIEW 01   05  4.00301   06  RCL 00   07  ST+ Y   08  *   09  STO M   10  3.002   11  LASTX   12  *   13  STO N   14  STO O   15  LBL 02   16  RCL IND O   17  XEQ IND X   18  LBL 03   19  STO IND M   20  DSE M   21  GTO 03   22  RCL 00   23  ENTER^   24  X^2   25  +   26  DSE X   27  ST+ M   28  DSE O   29  GTO 02 30  4   31  RCL 00   32  ST+ Y   33  *   34  STO M   35  RCL 00   36  LASTX   37  .1   38  %   39  +   40  +   41  STO O   42  LBL 04   43  RCL O   44  RCL 00   45  -   46  RCL IND X   47  X<> IND O   48  STO IND Y   49  LBL 05   50  RCL IND N   51  XEQ IND X   52  ST- IND M   53  DSE M   54  DSE N   55  GTO 05   56  RCL O   57  RCL 00   58  ST+ N 59  -   60  RCL IND X   61  X<> IND O   62  STO IND Y   63  DSE O   64  GTO 04   65  CLX   66  RCL 00   67  ST- M   68  X^2   69   E3   70  /   71  +   72   E3   73  /   74  RCL M   75  +   76  1   77  +   78  REGSWAP   79  RCL 00   80  RCL 00   81  LASTX   82  +   83  XEQ "LS3"   84  R^   85  INT   86  STO 00   87  STO M 88  STO N   89  ST+ N   90  STO O   91  ENTER^   92  X^2   93  +   94   E3   95  /   96  ISG X   97  RCL 00   98  3   99  * 100  ST+ O 101  + 102   E3 103  / 104  1 105  + 106  REGSWAP 107  X<>Y 108  X=0? 109  GTO 07 110  CLX  111  LBL 06 112  RCL IND M 113  ENTER^ 114  X<> IND N 115  - 116  RCL IND O 117  * 118  ST- IND M 119  ABS 120  + 121  DSE O 122  DSE N 123  DSE M 124  GTO 06 125   E-8 126  X M 142  CLA 143  CLD 144  END

( 241 bytes / SIZE n2+4n+1 )

 STACK INPUTS OUTPUTS X / Sum | fi |

-The same example gives - in 5mn37s -   | f1 | + | f2 | + | f3 | + | f4 | = 2 10-8  in X-register  and the solution in registers  R01 , R02 , R03 , R04

namely:   x1 =  4.266540475 , x2 =  1.353632241 ,  x3 = 3.548526764 ,  x4 = 6.831300511

3°) A short In/Out routine:

If flag F02 is clear, "IO" stores data.
If flag F02 is set, "IO" displays data.

RCL d  may be replaced by  RCLFLAG
STO d  may be replaced by  STO FLAG

-The append character is denoted  ~

 01  LBL "IO"  02  CF 22 03  CF 23 04  1 05  LBL 01       06  RCL d  07  FIX 0 08  CF 29 09  " "  10  ARCL Y 11  STO d 12  RDN 13  FS? 02 14  GTO 03 15  "~?"  16  PROMPT  17  FC?C 22 18  FS? 23 19  FS? 30 20  GTO 02       21  FS? 23 22  ENTER^ 23  FS?C 23 24  ASTO X 25  STO IND Z  26  CLX 27  SIGN 28  + 29  ISG Y 30  CLA 31  GTO 01       32  LBL 02 33  DSE X 34  X<>Y 35  - 36  CHS 37  DSE L 38  LASTX        39   E3 40  / 41  + 42  CLA 43  AOFF 44  RTN 45  LBL 03 46  "~="  47  ARCL IND Y 48  AVIEW 40  1 50  + 51  ISG Y 52  GTO 01 53  CLA 54  END

( 91 bytes )

 STACK CF02 INPUTS OUTPUTS X bbb bbb.eee
 STACK SF02 INPUTS OUTPUTS X bbb.eee /

Example:        Suppose you have to store  7 , 12 , 41 , "AAA" , "HP41"   into registers  R04 , R05 , ... etc ...

-Key in      CF 02

4    XEQ "IO"    the HP-41  displays     1?
7        R/S          the HP-41  displays     2?
12       R/S          the HP-41  displays     3?
41       R/S          the HP-41  displays     4?       switch to alpha mode and
"AAA"   R/S          the HP-41  displays     5?
"HP41"  R/S          the HP-41  displays     6?

simply press    R/S     and   you'll get   4.008     in X-register:  this is the control number  bbb.eee   of your data in registers R04 thru R08

-To obtain the control number needed for "LS", simply add  r.10-5 to the previous result

( r = the number of rows of the matrix = the number of equations of the system )

>>>  If you store a number like "PI" set flag F22 before pressing the R/S key

-If you want to view the contents of registers R04 to R08  press  SF 02    4.008    XEQ "IO"
and the HP-41 will display successively:

1 = 7
2 = 12
3 = 41
4 = AAA
5 = HP41

( if you set flag F21 , the calculator will stop at each  AVIEW )

>>>You can also use "IO" to store data into registers  R04 , R07 , R10 , ... etc ...

Key in   CF 02   4.00003  XEQ "IO"  and so on , but in this case, simply ignore the last prompt because the control number would be completely wrong!

-Likewise, to view registers   R04 , R07 , R10 , R13 , R16   for instance , key in   SF 02   4.01603  XEQ "IO"   ...

References:

  Francis Scheid - "Numerical Analysis" - ( McGraw-Hill )  ISBN:  07-055197-9
  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4