I'm trying to find the steady state temperature distribution in the infinite region outside a sphere of unit radius centred on the origin, where the temperature takes the value $u = f(\theta$) on the spherical surface, $((r ,\theta , \phi)$ being spherical polar coordinates), and $u \to 0$ as $r \to \infty$.
My textbook tells me that the solution is
$$u(r, \theta) = \sum\limits_{n = 0}^\infty \dfrac{A_n}{r^{n + 1}}P_n(\cos(\theta)),$$
where $$A_n = \dfrac{(2n + 1)}{2} \int_0^\pi f(\theta)P_n(\cos(\theta))\sin(\theta) \ d \theta,$$
where $P_n$ is the $n^{th}$ Legendre Polynomial, but I don't understand how it came to this solution.
I would greatly appreciate it if people could please take the time to explain how this problem is solved, so that I understand how it is done for similar types of problems.
The proposed solution is a classical expansion for the solution of Laplace equation with axial boundary conditions. In this symmetry, as noted in the comments, the eqation to be solved can be written as \begin{equation} \dfrac{\partial}{\partial{r}} \left( r^2 \dfrac{\partial{u}}{\partial{r}} \right) + \dfrac{1}{ \sin(\theta)} \dfrac{\partial}{\partial{\theta}}\left( \sin(\theta)\dfrac{\partial{u}}{\partial{\theta}} \right) = 0 \end{equation} Here is a brief description of the method. Solutions obtained by separation of the variables: $u(r,\theta)=R(r)\Theta(\theta)$ must verify \begin{equation} \frac{1}{R(r)}\dfrac{\partial}{\partial{r}} \left( r^2 \dfrac{\partial{R(r)}}{\partial{r}} \right)=-\frac{1}{\Theta(\theta)}\dfrac{1}{ \sin(\theta)} \dfrac{\partial}{\partial{\theta}}\left( \sin(\theta)\dfrac{\partial{\Theta(\theta)}}{\partial{\theta}} \right)=A \end{equation} where $A$ is a constant. The $\Theta$-ODE solutions are Legendre functions of $\cos\theta$. The non-divergence condition for $0\le \theta\le 2\pi$ implies \begin{equation} \frac{1}{2}\sqrt{1+4A}-\frac{1}{2}=\ell \end{equation} where $\ell$ is a positive integer, thus $\Theta(\theta)=b_\ell P_\ell (\cos\theta)$ (where $b_\ell$ is a constant). Then, the solution for $R$ which does not diverge for $r\to\infty$ is $R(r)=b_\ell r^{-\frac{1}{2}\sqrt{1+4A}-\frac{1}{2}}=b_\ell r^{-\ell-1}$. A general solution for the Laplace equation can thus be written as a superposition \begin{equation} u(r,\theta)=\sum_{\ell=0}a_lr^{-\ell-1}P_\ell \left( \cos\theta \right) \end{equation} Finally the constants $a_\ell$ are determined by the boundary condition $u(1,\theta)=f(\theta)$. Using the orthogonality properties of the Legendre polynomials, we obtain the given result.