Finding the memory polynomial model coefficients
Memory polynomial system identification is a linear regression problem.  Suppose you sent a $N$-element sequence $x(n)$ through a system and obtained an $N$-element output sequence $y(n)$.  If the system were described by a memory polynomial, then your model would be...

$y(n) = \sum_{q=0}^{Q}[x(n-q)a(q)+x(n-q)|x(n-q)|^2 b(q)] + \epsilon(n)$

Here $a(q)$ and $b(q)$ are coefficient of the models and $\epsilon(n)$ is an error term.  This is a linear model in $a$ and $b$, and can be expressed in matrix form as...

$\begin{bmatrix}y(Q)\\y(Q+1)\\...\\y(N-1)\end{bmatrix} = \begin{bmatrix}x(Q)&amp;...&amp;x(0)&amp;x(Q)|x(Q)|^2&amp;...&amp;x(0)|x(0)|^2\\x(Q+1)&amp;...&amp;x(1)&amp;x(Q+1)|x(Q+1)|^2&amp;\ldots&amp;x(1)|x(1)|^2\\\ldots&amp;\ldots&amp;\ldots&amp;\ldots&amp;\ldots&amp;\ldots\\x(N-1)&amp;\ldots&amp;x(N-1-Q)&amp;x(N-1)|x(N-1)|^2&amp;\ldots&amp;x(N-1-Q)|x(N-1-Q)|^2\\\end{bmatrix}\begin{bmatrix}a(0)\\...\\a(Q)\\b(0)\\...\\b(Q)\end{bmatrix}+\begin{bmatrix}\epsilon(0)\\\epsilon(1)\\...\\\epsilon(N-1)\end{bmatrix}$

You can solve for $a$ and $b$ using any of the linear least squares functions in MATLAB, Python, etc.