def solve_leq(A,B) { N = size(A)[0]; for ( J = 0; J < N; J++ ) { for ( I = J; I < N; I++ ) if ( A[I][J] != 0 ) break; if ( I == N ) return 0; if ( I != J ) { for ( L = 0; L < N; L++ ) { T = A[J][L]; A[J][L] = A[I][L]; A[I][L] = T; } T = B[J]; B[J] = B[I]; B[I] = T; } for ( K = J+1; K < N; K++ ) { M = A[K][J]/A[J][J]; for ( L = J; L < N; L++ ) A[K][L] -= M*A[J][L]; B[K] -= M*B[J]; } } for ( J = N-1; J >= 0; J-- ) { for ( K = J+1, S = 0; K < N; K++ ) S += A[J][K]*B[K]; B[J] = (B[J]-S)/A[J][J]; } return 1; } A=newmat(3,3,[[1,2,3],[4,5,6],[7,8,8]]); B=newvect(3,[1,2,3]); printf("ans is %a\n",matrix_inverse(A)*B)$ solve_leq(A,B)$ print(B)$ end$