Solution (source code)

= Solution

Set $\Phi=R(r)P(\theta)$ and multiply the <Laplace equation> by $r^2/(RP)$. <Separation of variables> gives
$$
\frac{(r^2R')'}R=-\frac{(\sin\theta P')'}{\sin\theta P}=\lambda.
$$
With $\xi=\cos\theta$, the angular equation is $[(1-\xi^2)P_\xi]' +\lambda P=0$. Regularity at both poles selects the <Legendre polynomials> $P_n$ and $\lambda=n(n+1)$. The radial <Cauchy-Euler equation> has solutions $r^n,r^{-n-1}$. Expanding the angular data in this complete orthogonal system therefore gives
$$
\Phi(r,\theta)=\sum_{n=0}^\infty(A_nr^n+B_nr^{-n-1})P_n(\cos\theta),
$$
with the usual convergence required of a <harmonic> expansion on the region considered. Angular Legendre functions singular at the poles are excluded.

For $0<a<b$, the outer data have expansion $b^2\cos^2\theta=b^2/3+(2b^2/3)P_2(\cos\theta)$, while the inner data are pure $P_0$. Thus only degrees zero and two are needed. A radial interpolation matching each pair of boundary values gives
$$
\boxed{\Phi(r,\theta)=\frac{a^3(b-r)+(b^3/3)(r-a)}{r(b-a)}
+\frac{2b^2}{3}\frac{r^2-a^5/r^3}{b^2-a^5/b^3}P_2(\cos\theta).}
$$
The first term is a <linear combination> of $1,r^{-1}$ and takes values $a^2,b^2/3$ at the two radii. The second is a combination of $r^2,r^{-3}$ and takes values zero and $2b^2/3$. These checks establish both boundary conditions; each term is <harmonic>. The <maximum principle for harmonic functions> gives uniqueness. This is an example of <axisymmetric harmonic interpolation between concentric spheres>.