 Solution

For an axisymmetric harmonic function, put $\Psi =R(r)\Theta (\theta )$ in the Laplacian in spherical coordinates and multiply by $r^2/(R\Theta )$. Separation of variables givesWith $u=\cos \theta$, the angular equation becomes the Legendre differential equation. Regularity at both poles selects the Legendre polynomials $P_{\ell }(u)$ with integer $\ell \ge 0$; the singular second solutions are excluded. The radial Euler-Cauchy equation has solutions $r^{\ell }$ and $r^{-\ell -1}$. Superposition therefore giveswith coefficients chosen to ensure convergence and the required boundary regularity. The axisymmetric harmonic extension from a sphere gives these coefficients directly from boundary data.