The GP equation is
$$i\frac{\partial u}{\partial z}+\nabla^2u+|u|^2u+\int e^{-[(x-x')^2+(y-y')^2]}|u(x',y')|^2 dx'dy'u(x,y)=0$$
with Neumann boundary condition. The initial condition is given by a Gaussian function, and $x, y\in [-10,10], z>0$.
Question: How do we deal with partial differential equations with a convolution integral?