docs: public double[][] solve(double[][] B) {
          if (B.length != 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[0].length;
          double[][] X = Matrix.Submatrix(B, piv, 0, nx - 1);
  
          // 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 X;
      }