48$. The difference between the values used for $N_x$ and $N_y$ follows from the observation that if the solution in the $y$-direction is smooth enough, then only few low wave-vectors are exited. The value of the parameter $y_{\rm max}$ fixes the $y$-range where the solution is assumed different form zero, since in the numerical algorithm is implicitly assumed that $$P(x,y)\equiv m(x,y) = 0 \qquad |y|> y_{\rm max}.$$ This explain the rather large value used. The number of iterations necessary to reach a mean square error on $q$, $P$ and $m$ of order $O(10^{-6})$ depends on the initial guess of $q(x)$ but it is typically of few hundreds.

In Figure \[fig:qx-pdq.T0.6-1.000\] are shown the order parameter $q(x)$ and the overlap probability distribution function $P(q)$ at $T=0.6$ computed for increasing $x$-resolution and $x_0=1$. As expected the agreement between the numerical and the perturbative solutions increases with the number of $N_x$ of $x$-grid points. However, the convergence is not uniform: it is rather fast far from $x_{\rm max}$ and much slower for $x\simeq x_{\rm max}$, see the inset of Figure \[fig:qx-pdq.T0.6-1.000\] panel (a). This is not unexpected because for $x=x_{\rm max}$ the derivative of the order parameter $q^{(1)}(x)$ has a cusp: $$\label{eq:cusp}
\lim_{x\to 