def METHOD_NAME(N, w, rows):
    if N == 1:
        subset = w * np.ones([rows, 1])
    else:
        if w == 0:
            subset = np.zeros([rows, N])
        elif w == 1:
            subset = np.eye(N)
        else:
            # initialize submatrix
            subset = np.empty([rows, N])
            # starting row of submatrix
            row_start = 0
            # iterate by polynomial order and fill the multiindex submatrices
            for k in range(0, w + 1):
                # number of rows of the submatrix
                sub_rows = PolynomialBasis._setsize(N - 1, w - k)
                # update until row r2
                row_end = row_start + sub_rows - 1
                # first column
                subset[row_start:row_end + 1, 0] = k * np.ones(sub_rows)
                # subset update --> recursive call
                subset[row_start:row_end + 1, 1:] = \
                    PolynomialBasis.METHOD_NAME(N - 1, w - k, sub_rows)
                # update row indices
                row_start = row_end + 1
    return subset