Solve Laplace's equation in spherical coordinates, $\nabla^2 u(r,\theta,\phi)=0$, in the general case. To find the equilibrium temperature distribution of a sphere I am required to solve the following PDE BVP (Laplace's equation in spherical coordinates).

$$\nabla^2 u(r,\theta,\phi)=0, 0<r<a,0<\theta<2\pi,0<\phi<\pi$$
  $$u(a,\theta,\phi)=f(\theta,\phi), 0<\phi<\pi$$

Note that I am required to solve this problem for the general case.
Work thus far
Note that in spherical coordinates $\nabla^2 u(r,\theta,\phi)=u_{rr}+\frac{2}{r}u_r+\frac{1}{r^2}[u_{\phi \phi}+\cot(\phi)u_{\phi}+\csc^2(\phi)u_{\theta \theta}]=0$
We use separation of variables letting $u(r,\theta,\phi)=w(r)v(\theta,\phi)$.
Substituting in
\begin{align}
0&=u_{rr}+\frac{2}{r}u_r+\frac{1}{r^2}[u_{\phi \phi}+\cot(\phi)u_{\phi}+\csc^2(\phi)u_{\theta \theta}] \\
0&=w''v+\frac{2}{r}w'v+\frac{1}{r^2}[wv_{\phi \phi}+\cot(\phi)wv_{\phi}+\csc^2(\phi)wv_{\theta \theta}] \\
\frac{r^2}{w}(-w''-\frac{2}{r}w')&=\frac{1}{v}[v_{\phi \phi}+\cot(\phi)v_{\phi}+\csc^2(\phi)v_{\theta \theta}]=-\lambda \\
\end{align}
to obtain the following ODEs
\begin{align}
0&=r^2w''+2rw'-\lambda w\\
0&=v_{\phi \phi}+\cot(\phi)v_{\phi}+\csc^2(\phi)v_{\theta \theta}+\lambda v\\
\end{align}
We now do another separation of variables $v(\theta,\phi)=g(\theta)h(\phi)$. We get the following separation
$$\frac{h''+\cot(\phi)h'+\lambda h}{\csc^2(\phi)h}=\frac{-g''}{g}=\mu $$
and obtain the following ODEs
\begin{align}
0&=g''+\mu g\\
0&=h'' + \cot(\phi)h'+\lambda h-\mu \csc^2(\phi)h\\
\end{align}
Solving the eigenvalue problem $g''+\mu g=0$ with periodic boundary conditions $g(0)=g(2\pi)$ and $g'(0)=g'(2\pi)$, we obtain the eigenvalues $\mu_n=n^2$ and the eigenvectors $g_n(\theta)=\cos(n\theta)$ or $\sin(n\theta)$ for $n=0,1,2,3,...$.
What I am confused about
What order do I solve the eigenvalue problems to find the eigenvalues $\mu,\lambda$?
How do I transform the ODE $0=h'' + \cot(\phi)h'+\lambda h-\mu \csc^2(\phi)h$ into a easier ODE to solve.
Note that following posts are not relevant as 


*

*Radial symmetry is used

*Derive the formula for $\nabla^2 u(r,\theta,\phi)$
 A: Assume
$$
  x=r\sin\theta\cos\phi,\;\;\; y=r\sin\theta\sin\phi,\;\;\; z=r\cos\theta, \\
   0 \le r < \infty, \;\;\; 0 \le \theta \le \phi,\;\;\; 0 \le \phi \le 2\pi.
$$
Because the spherical coordinate system is an orthogonal coordinate system (meaning that the coordinate curves where only one variable varies are mutually orthogonal where they meet,) the Laplacian is determined by the metric scale factors
$$
          m_r = 1,\;\; m_{\theta}=r,\;\; m_{\phi} = r\sin\theta.
$$
The Laplacian is
\begin{eqnarray*}
 & & \nabla^{2}F(r,\phi,\theta) \\
   & = & \frac{1}{m_{r}m_{\phi}m_{\theta}}\left[
      \frac{\partial}{\partial r}\left(\frac{m_{\phi}m_{\theta}}{m_{r}}
            \frac{\partial F}{\partial r}\right)
     +\frac{\partial}{\partial\phi}\left(\frac{m_{r}m_{\theta}}{m_{\phi}}
            \frac{\partial F}{\partial\phi}\right)
     +\frac{\partial}{\partial\theta}\left(\frac{m_{r}m_{\phi}}{m_{\theta}}
            \frac{\partial F}{\partial\theta}\right)
       \right]
       \\
%   & = & \frac{1}{r^{2}\sin\theta}\left[
%      \frac{\partial}{\partial r}\left(r^{2}\sin\theta
%            \frac{\partial F}{\partial r}\right)
%     +\frac{\partial}{\partial\phi}\left(\frac{1}{\sin\theta}
%            \frac{\partial F}{\partial \phi}\right)
%     +\frac{\partial}{\partial\theta}\left(\sin\theta
%            \frac{\partial F}{\partial\theta}\right)
%       \right]
%       \\
   & = & \frac{1}{r^{2}}\frac{\partial}{\partial r}
        \left(r^{2}\frac{\partial F}{\partial r}\right)
     + \frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}F}{\partial\phi^{2}}
     + \frac{1}{r^{2}\sin\theta}
        \frac{\partial}{\partial\theta}
        \left(\sin\theta\frac{\partial F}{\partial\theta}\right).
\end{eqnarray*}
Standard separation of variables equations are obtained by setting $F(r,\theta,\phi)=R(r)\Theta(\theta)\Phi(\phi)$. The separated equations involve a separation parameter $m=0,\pm 1,\pm 2,\cdots$ imposed by periodicity in $\phi$, and a separation parameter $l$ determined by the ODE in $\theta$. The separated equations are:
\begin{eqnarray*}
       -\frac{d^2}{d\phi^2}\Phi & = & m^2\Phi, \\
   -\frac{d}{d\theta}\left(\sin\theta\frac{d\Theta}{d\theta}\right)
       +\frac{m^2}{\sin\theta}\Theta & = & l\sin\theta \Theta, \\
   -\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+l R & = & 0.
\end{eqnarray*}
The equation in $\theta$ is the Associated Legendre equation, and this equation determines the possible values of $l$. The radial equation is determined by the requirement that $R$ be bounded near $r=0$. The equation in $\theta$ is transformed to standard Legendre form by setting $x=-\cos\theta$ so that $x$ varies between $-1$ and $1$ as $\theta$ varies between $0$ and $\pi$. Then
$X$ defined by $X(-\cos\theta)=\Theta(\theta)$ satisfies the Associated Legendre Equation:
$$
       -\frac{d}{dx}(1-x^2)\frac{dX}{dx}+\frac{m^2}{1-x^2}X=l X.
$$
For a given $m$, the values $l$ for which there are bounded solutions $X$ are those given by $l=n(n+1)$ where $n=m,m+1,m+2,\cdots$. The resulting bounded solutions are scalar multiples of
$$
           e^{im\phi}r^{n}\frac{d^{n+m}}{dx^{n+m}}(1-x^2)^{n}.
$$
