Solution (source code)

= Solution

Write the <scalar field> as $u=R(r)Y(\theta,\phi)$, where $R=r+r^{-2}$ and $Y=\sin\theta\cos\phi$. Its radial contribution to the <Laplacian> in <spherical coordinates> is
$$
\frac{Y}{r^2}(r^2R')'=Y\left(\frac2r+\frac2{r^4}\right)=\frac{2R}{r^2}Y.
$$
For the angular contribution, $Y_\theta=\cos\theta\cos\phi$ and $Y_{\phi\phi}=-\sin\theta\cos\phi$, so
$$
\frac1{\sin\theta}\partial_\theta(\sin\theta\,Y_\theta)
+\frac{Y_{\phi\phi}}{\sin^2\theta}
=\frac{\cos\phi}{\sin\theta}(\cos^2\theta-\sin^2\theta-1)=-2Y.
$$
Combining the radial and angular terms proves
$$
\boxed{\nabla^2u=\frac{2R}{r^2}Y-\frac{2R}{r^2}Y=0\qquad(r>0).}
$$
The calculation uses $0<\theta<\pi$ to avoid coordinate singularities. In <Cartesian coordinates>, $u=x(1+r^{-3})$, which extends smoothly across the polar axis for $r>0$. It therefore satisfies <Laplace's equation> there as well. The origin is a genuine singularity and is excluded.