Use of Legendre's equation. For some weeks have been studying Legendre polynomial as a solution to this equation. 
$$ (1-x^2)\frac{d^2}{dx^2}f(x)-2x\frac{d}{dx}f(x)+n(n+1)f(x)=0.$$
I've found them very interesting to learn from purely mathematical perspective but I haven't come across their specific use.  Why are they so important and any example suggesting the use of Legendre polynomial would be appreciated.
Thank you.
 A: One of the greatest applications not only of the Legendre polynomials, but of orthogonal polynomials in general, is Quantum Mechanics. A great example is the one electron atom:
The one electron atom is a three dymensional system containing two particles, a nucleus and an electron, moving under influence of a mutual Coulomb attraction and bounded togheter by that attraction. Let the mass of the nucleus be $m_1$ and its charge be $+Ze$ ($Z=1$ for neutral hydrogen and $Z=2$ for singly ionized helium, etc.). Designate the mass of the electron by $m_2$ and its charge by $-e$. Then the Coulomb potential energy of the system is
$$
V(x_1,y_1,z_1,x_2,y_2,z_2) = - \frac{Z e^2}{\sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2 + (z_2 - z_1)^2}}
$$
where $x_1,y_1,z_1$ are the rectangular coordinates of the nucleus and $x_2,y_2,z_2$ are the rectangular coordenates of the electron. The time independent Schrödinger equation for this case is
$$
-\frac{\hbar^2}{2 m_1} \Delta_1 \psi_T - \frac{\hbar^2}{2 m_2} \Delta_2 \psi_T + V(x_1,\ldots,z_2) = E_T \psi_T,
$$
where $\Delta_i = \frac{\partial^2}{\partial x_i^2} + \frac{\partial^2}{\partial y_i^2} + \frac{\partial^2}{\partial z_i^2}$ is the Laplacian operator.
Taking the change of variables
\begin{align}
x &= \frac{m_1 x_1 + m_2 x_2}{m_1 + m_2}\\
y &= \frac{m_1 y_1 + m_2 y_2}{m_1 + m_2}\\
z &= \frac{m_1 z_1 + m_2 z_2}{m_1 + m_2}\\
\end{align}
the rectangular coordinates of the center of mass of the system, and
\begin{align}
r \sin \theta \cos \phi &= x_2 - x_1\\
r \sin \theta \sin \phi &= y_2 - y_1\\
r \cos \theta &= z_2 - z_1
\end{align}
the spherical coordenates of the electron relative to the nucleus, the Schrödinger equation can be transformed to
$$
-\frac{\hbar^2}{2 M} \Delta_{(x,y,z)} \psi_T - \frac{\hbar^2}{2 \mu} \Delta_{(r, \theta, \phi)} \psi + V(r) \psi_T = E_T \phi_T
$$
where
$$
M = m_1 + m_2, \quad \mu = \frac{m_1 m_2}{m_1 + m_2}, \quad V(r) = -\frac{Z e^2}{r},
$$
$\Delta_{(x,y,z)}$ is the Laplace operator in the cartesian coordinates $(x,y,z)$ and $\Delta_{(r,\theta,\pi)}$ is the Laplace operator in the spherical coordinates $(r, \theta, \phi)$.
The separation of this time independent Schrödinger equation is effected by looking for solutions of the form
$$
\psi_T(x,y,z,r,\theta,\phi) = \psi_{CM}(x,y,z) \psi(r,\theta,\psi).
$$
Substituting this into the equation and dividing by $\psi_T = \psi_{CM} \psi$ gives
$$
\left[-\frac{\hbar^2}{2 M} \frac{\Delta_{(x,y,z)} \psi_{CM}}{\psi_{CM}}\right] + \left[-\frac{\hbar^2}{2 \mu} \frac{\Delta_{(r,\theta,\psi)} \psi}{\psi} + V(r)\right] = E_T.
$$
The first term on the left is a function only of $x,y,z$ and the second term is a function only of $r,\theta,\phi$. Since the sum of the two terms is equal to a constant $E_T$, each term must separately be equal to a constant. Thus we obtain two equations,
$$
\left[-\frac{\hbar^2}{2 M} \frac{\Delta_{(x,y,z)} \psi_{CM}}{\psi_{CM}}\right] = E_{CM}
$$
and
$$
\left[-\frac{\hbar^2}{2 \mu} \frac{\Delta_{(r,\theta,\psi)} \psi}{\psi} + V(r)\right] = E
$$
where $E_{CM} + E = E_T$.
The first equation is the time independent Schrödinger equation for a particle of mass $M$, whose coordinates are those of the center of mass of the atom, and which is moving with total energy $E_{CM}$ in a region of zero potential. The equation describes the translational motion of the atom. The second equation is the time independent Schrödinger equation describing the motion of the electron relative to the nucleus. Since $E = E_T - E_{CM}$, where $E_T$ is the total energy of the atom and $E_{CM}$ is its energy of translational motion, $E$ is the total energy of relative motion -that is, the total energy of the atom in a frame of reference in which the center of mass is stationary-. This equation is formally identical with the time independent Schrödinger equation for an electron of mass $\mu$, moving with coordinates $r,\theta,\phi$ about an infinitley massive nucleus that remains stationary.
Given the form of the Coulomb potential and, in fact, whenever the potential $V$ is any function of the radial coordinate $r$ only, it is posible to separate the equation into three ordinary differential equations. To do this, we look for solutions of the form
$$
\psi(r,\theta,\phi) = R(r) \Theta(\theta) \Phi(\phi)
$$
Inserting this into the second equations and dividing through by $-\frac{R\Theta\Phi \hbar^2}{2\mu}$, we can separate variables and obtain the three equations
$$
\frac{d^2 \Phi}{d \phi^2} = -m^2 \Phi \tag{Harmonic Oscillator}
$$
$$
-\frac{1}{\sin \theta} \frac{d}{d\theta} \left(\sin \theta \frac{d \Theta}{d \theta}\right) + \frac{m^2 \Theta}{\sin^2 \theta} = \alpha \Theta
$$
and
$$
\frac{1}{r^2} \frac{d}{dr}\left(r^2 \frac{d R}{d r}\right) + \frac{2 \mu}{\hbar^2}[E - V(r)]R = \alpha \frac{R}{r^2}
$$
where $m$ and $\alpha$ are the separation constants.
The solution for the azimuthal equation is easily verified to be $\Phi(\phi) = e^{i m \phi}$, and it must be single-valued. Thus we must make sure that $\Phi(0) = \Phi(2\pi)$, or $|m| = 0,1,2,3,\ldots$, and the set of solutions which are acceptable are
$$
\Phi_m(\phi) = e^{i m \phi}
$$
where $m \in \mathbb{Z}$.
For the polar equation, it's convenient to change variables from $\theta$ to $\xi$, where $\xi = \cos \theta$. Then the equation reads
\begin{align}
\frac{d}{d \xi}\left\{(1-\xi^2)\frac{d \Theta(\xi)}{d \xi}\right\} + \left\{\alpha - \frac{m^2}{1-\xi^2}\right\}\Theta(\xi) = 0 \\ {}\tag{Associated Legendre Equation}
\end{align}
which, considering the behavior at $\xi = \pm 1$, can be transformed using the change of variables
$$
\Theta(\xi) = (1 - \xi^2)^{|m|/2} P(\xi)
$$
to the differential equation
$$
\frac{d}{d\xi}\left[(1-\xi^2) \frac{d P}{d\xi}\right] + \alpha P(\xi) = 0 \tag{Legendre equation}
$$
The solutions of this equations are the Legendre functions, and the condition that $\Theta(\xi)$ is everywhere finite implies that $P$ must be a polynomial. This happens only when $\alpha$ has the value
$$
\alpha = l(l+1)
$$
where $l =|m|,\,|m|+1,\,|m|+2,\ldots$. The polynomials $P_l(\xi)$ are the Legendre polynomials, and the functions $\Theta_{lm}(\xi) = (1-\xi^2)^{|m|/2} P_l(\xi)$ are the associated Legendre polynomials.
Lastly, taking the new variable $\rho = 2 \beta r$ where $\beta^2 = - \frac{2\mu E}{\hbar^2}$ and using $\gamma = \frac{\mu Z e^2}{\hbar^2 \beta}$, the radial equation becomes
$$
\frac{1}{\rho^2} \frac{d}{d \rho} \left(\rho^2 \frac{d R(\rho)}{d \rho}\right) + \left\{-\frac{1}{4} - \frac{l(l+1)}{\rho^2} + \frac{\gamma}{\rho}\right\}R(\rho) = 0
$$
Considering the behavior of the equation for $\rho \to 0$ and $\rho \to \infty$, we use the transformation
$$
R(\rho) = e^{-\rho/2} \rho^l L(\rho)
$$
and reduce the equation to
$$
\rho \frac{d^2 L}{d \rho^2} + (1-x)\frac{d L}{d \rho} + \gamma L = 0 \tag{Laguerre Equation}
$$
To prevent $R(\rho)$ from diverging at $\rho \to \infty$, it is found that $\gamma = n$, where $n$ is one of the integers $n = l+1,\,l+2,\ldots$, and the acceptable radial functions are $R_{nl}(\rho) = e^{-\rho/2}\rho^l L_n(\rho)$. The polynomials $L_n(\rho)$ are known as Laguerre polynomials and the functions $R_{nl}$ are the associated Laguerre functions.
In a nutshell, the one electron atom has as eigenfunctions
$$
\psi_{nlm}(r,\theta,\phi) = e^{-\beta r}(2\beta r)^l \sin^{|m|}\theta \,L_n(2 \beta r) P_l (\cos \theta) e^{i m \phi},
$$
quantum numbers
\begin{align}
n &= 1,2,3,\ldots\\
l &= 0,1,2,\ldots, n-1\\
m &= -l,-l+1,\ldots,0,\ldots,l-1,1
\end{align}
and allowed values for the toatal energy
$$
E_n = - \frac{\mu Z^2 e^4}{2\hbar^2 n^2}, \quad n=1,2,3,4,\ldots
$$
which are the energy levels for the single bound electron atom.
We can conclude that quantization is a mathematical consequence of the Schrödinger Equation, and is realized by means of orthogonal polynomials!
A: Roots of such polynomials are used through Gaussian quadrature, which aims at computing integrals accurately.
For instance the second Legendre polynomial, namely
$$ P_2(X) = \frac{1}{2}\left(3X^2-1\right) $$
has roots $\pm \frac{1}{\sqrt{3}}$. Gaussian quadrature tells you that 
$$f\left(\frac{1}{\sqrt{3}}\right) + f\left(-\frac{1}{\sqrt{3}}\right)$$
is a good approximation for
$$I = \int_{-1}^1 f(x) dx$$
In fact it is exact for polynomials up to the 3th order !
You can get more info there as well.
