Solution (source code)

= Solution

The first-order condition $Q=0$, combined with $b_0=-i\psi_{0,y}$, gives
$$
b_1=-i\psi_{1,y}-\frac i\omega\psi_{0,x}+\frac1\omega\psi_{0,yyy}.
$$
Differentiate this expression with respect to $y$ and insert it into $S=0$. The $\psi_1$ terms cancel, leaving
$$
\begin{aligned}
0&=-i\omega\psi_{1,yy}-\omega\left[-i\psi_{1,yy}-\frac i\omega\psi_{0,xy}+\frac1\omega\psi_{0,yyyy}\right]+i\psi_{0,xy}\\
&=2i\psi_{0,xy}-\psi_{0,yyyy}.
\end{aligned}
$$
Thus the compatibility <partial differential equation> is
$$
\boxed{2\psi_{0,xy}+i\psi_{0,yyyy}=0.}
$$
It is exactly $\partial_y(2\psi_{0,x}+i\psi_{0,yyy})=0$. Integrating with respect to $y$ therefore gives
$$
\boxed{2\psi_{0,x}+i\psi_{0,yyy}=f(x),}
$$
where the integration constant may be any function of $x$. Conversely, differentiating this last equation with respect to $y$ recovers the compatibility equation.