docs: public Matrix solve (Matrix B) {
        if (B.getRowDimension() != m) {
           throw new IllegalArgumentException("Matrix row dimensions must agree.");
        }
        if (!this.isNonsingular()) {
           throw new RuntimeException("Matrix is singular.");
        }
  
        // Copy right hand side with pivoting
        int nx = B.getColumnDimension();
        Matrix Xmat = B.getMatrix(piv,0,nx-1);
        double[][] X = Xmat.getArray();
  
        // Solve L*Y = B(piv,:)
        for (int k = 0; k < n; k++) {
           for (int i = k+1; i < n; i++) {
              for (int j = 0; j < nx; j++) {
                 X[i][j] -= X[k][j]*LU[i][k];
              }
           }
        }
        // Solve U*X = Y;
        for (int k = n-1; k >= 0; k--) {
           for (int j = 0; j < nx; j++) {
              X[k][j] /= LU[k][k];
           }
           for (int i = 0; i < k; i++) {
              for (int j = 0; j < nx; j++) {
                 X[i][j] -= X[k][j]*LU[i][k];
              }
           }
        }
        return Xmat;
     }