Solving lin.eq and matrix inversions is limited to matrices with sizes up to 8x8, by J.Fossy Weinzinger The hp15c was a highlight of its time. It is noteworthy that the original nut code of hp15c is apparently fit for an extension of its RAM. As far as I can say it work best for all but the inverse of a matrix larger than 8x8. The largest matrix the original hp15c could store was 8x8 - so the algorithm was optimized for that maximum in size.

The reason why the original algorithm fail for larger matrices is that it has to reorder the rows of the matrix under certain circumstances. How the rows are reordered is stored in the signs of the diagonal elements of the matrix. In a register (7 bytes - 14 nibbles) the sign is stored in a nibble as 0 for + or 9 for -. A nibble hold 4 bits - so in the sign is place for additional 3 bits. These remaining 3 bits are used to store the original row offset. This only work up to 8 rows. If a matrix has more than 8 rows this will not work any more. :-(

The intermediate step used for matrix division, matrix inverse and determinant of a matrix is called LU decomposition. It is described in the "HP15C advanced functions handbook" - Section 4 - "Using Matrix Operations" - "Understanding the LU Decomposition" on page 96/97 of the original manual from 1982 or an page 82/83 of the reproduction from 2012.

Do not trust the results of matrix division, matrix inverse or determinant of a matrices larger than 8x8.

If you are interested in a easy counter-example:

---------------------------------------------
A(i,j) = 0 if j < (n - i + 1);
Otherwise A(i,j) = 1

b(i) = i
1 <= i,j <= n

A*x = b

=>

x(i) == 1 for all 1 <= i <= n

---------------------------------------------
A * x = b

0 0 0 0 ... 0 1 x(1) = 1
0 0 0 0 ... 1 1 x(2) = 2
...
0 1 1 1 ... 1 1 x(n-1) = n-1
1 1 1 1 ... 1 1 x(n) = n

---------------------------------------------

The following program, also available as a dump, will fill pre-dimensioned matrix A, set size of b, fill b and set result matrix to C.

f LBL D 001-42,21,14
RCL DIM A 002-45,23,11
STO 2 003-   44  2
g TEST 6 (x≠y?) 004-43,30, 6
g RTN 005-   43 32
f MATRIX 1 (R0,R1 := 0)006-42,16, 1
f LBL 0 007-42,21, 0
RCL 2 008-   45  2
RCL 0 009-   45  0
- 010-      30
1 011-       1
+ 012-      40
RCL 1 013-   45  1
g TEST 8 (x<y?) 014-43,30, 8
g CLx 015-   43 45
g TEST 1 (x>0?) 016-43,30, 1
1 017-       1
USER STO A (USER) 018u   44 11
GTO 0 019-   22  0
RCL 2 020-   45  2
1 021-       1
f DIM B 022-42,23,12
f MATRIX 1 (R0,R1 := 0)023-42,16, 1
f LBL 1 024-42,21, 1
RCL 0 025-   45  0
USER STO B (USER) 026u   44 12
GTO 1 027-   22  1
f RESULT C 028-42,26,13
g RTN 029-   43 32

Then use it like this:

<n> ENTER f DIM A
f D
RCL MATRIX b
RCL MATRIX A
/
f USER
f MATRIX 1
RCL C
RCL C
...
f USER


Do the test:
If 0<n<9 all elements of C will be equal to 1, but if n>8 the algorithm will fail and C will hold totally wrong values. Due to an overflow, the display will blink.