Solution (source code)

= Solution

For a separated <axisymmetric harmonic function> $\Psi=R(r)\Theta(\theta)$, multiplying the <Laplace equation> by $r^2/(R\Theta)$ gives
$$
\frac{(r^2R')'}{R}
=-\frac{(\sin\theta\,\Theta')'}{\sin\theta\,\Theta}=\lambda.
$$
The two sides depend on different variables, hence both equal a constant. Put $u=\cos\theta$ in the angular equation. Since $\sin\theta\,\Theta'=-\sin^2\theta\,d\Theta/du$, it becomes the <Legendre differential equation>
$$
\frac{d}{du}\left[(1-u^2)\frac{d\Theta}{du}\right]+\lambda\Theta=0.
$$
Regularity on both parts of the polar axis restricts $\lambda=\ell(\ell+1)$, with $\ell=0,1,\ldots$, and selects the <Legendre polynomial> $P_\ell(u)$. The other angular solution is singular at an endpoint.

The radial <Euler differential equation> is
$$
r^2R''+2rR'-\ell(\ell+1)R=0.
$$
Trying $R=r^s$ gives $s(s+1)=\ell(\ell+1)$, so $s=\ell$ or $s=-\ell-1$. <Superposition> of these separated <harmonic functions> gives
$$
\boxed{\Psi(r,\theta)=\sum_{\ell=0}^\infty
\left(A_\ell r^\ell+B_\ell r^{-\ell-1}\right)P_\ell(\cos\theta)}.
$$
The <Legendre polynomials> form the angular expansion for regular data with <rotational symmetry> about the polar axis; the radial domain and the <boundary conditions> determine which coefficients survive and where the expansion converges.