Solution
= Solution
Use the large-$z_0$ expansion
$$
\frac1{\sqrt{z_0^2+R^2}}
=\frac1{z_0}-\frac{R^2}{2z_0^3}+O(z_0^{-5}).
$$
The first term vanishes by the compatibility condition. For the given $p$, symmetry leaves
$$
\int p(x,y)\bigl((x-x_0)^2+(y-y_0)^2\bigr),dx,dy
=-2x_0\left(\int_{-\pi/2}^{\pi/2}x\sin x\,dx\right)\pi
=-4\pi x_0.
$$
Consequently
$$
\boxed{
u(x_0,y_0,z_0)
=-\frac{x_0}{z_0^3}+O(z_0^{-5})}.
$$