I am trying to reproduce the numerical solution of the so-called Parisi equations for the Sherrington-Kirkpatrick (SK) model. There are at least two main methods to solve these equations, and unfortunately I am confused by both.

First, the equations themselves. They are:

$$ q(x) = \int_{-\infty}^{\infty} dy \, P(x,y) m(x,y)^2 \,,$$

$$ \dot{m}(x,y) = - \frac{\dot{q}(x)}{2} \left[ m''(x,y) + 2 \beta x m(x,y) m'(x,y) \right] \,,$$

$$ \dot{P}(x,y) = \frac{\dot{q}(x)}{2} \left[ P''(x,y) - 2 \beta x \left[ m(x,y) P(x,y) \right]' \right] \,,$$

and the initial conditions are:

$$ m(1,y) = \tanh(\beta y), \qquad P(0,y) = \delta(y) \,.$$

Primes denote $y$ derivatives and a dot denotes an $x$ derivative. The domain is $x \in [0,1]$, $y \in [-\infty, \infty]$. In any numerical attempt, the range of $y$ is approximated as finite, $y \in [-y_{\text{max}}, y_{\text{max}}]$.

In each of the ways to solve these equations that I know about, one proceeds by first supplying a guess for $q(x)$, then solving the $m$ equation, then the $P$ equation. Then $q$ is updated according to the first equation. These steps are iterated until convergence, which typically takes hundreds of cycles .

The two main methods I know of are an integral kernel method by Nemoto and a pseudo-spectral method by Crisanti and Rizzo. In the Nemoto method kernels are used, and these require an initial guess for not only $q(x)$ but also $m(x,y), P(x,y)$. I have no idea what a good guess would be, and I have tried many and failed to get convergence. Presumably the guess should agree with the initial conditions for $m, P$ for $x=1, 0$ respectively. The Crisanti and Rizzo method takes a Fourier transform in the $y$-direction, and then numerically integrates in the $x$-direction using an Adams-Bashforth scheme. Implicit in their approach is that they are setting $P, m$ to have Dirichlet boundary conditions in the $y$-direction, although this is inconsistent with the initial condition for $m$. I have tried to reproduce their method using using both Mathematica's NDSolve, and by adapting this pseudo-spectral code by Trefethen for the KdV equation, but I cannot get either code to converge. From the shape of the curves it is clear I am not handling the boundary conditions properly.

I would greatly appreciate any guidance, as well as any references that might shed light or offer a more easy-to-understand numerical method to solve these equations.