Solution (source code)

= Solution

For a <monomial> trial $u_p=Ax^my^n$, the operator gives
$$
\mathcal Lu_p=2Amx^{m-1}y^{n+1}-Anx^{m+1}y^{n-1}.
$$
Choosing $m=n=2$ matches the two powers in the forcing, and $A=1/2$ matches their coefficients. A <particular solution> is therefore $u_p=x^2y^2/2$. Adding the homogeneous solutions gives
$$
\boxed{u(x,y)=\frac12x^2y^2+F(x^2+2y^2)}.
$$
The arbitrary <function> labels the values on the <characteristic curves>; it must have sufficient regularity for the desired classical solution.