Well, the best way I can think of to introduce a real physical problem where Laguerre ànd Legendre (even of the associated form!) differential equations play a role, is the (exact, that is, without taking spin into consideration) solution of the combined angular momentum / magnetic moment / energy eigenvalue-problem for the hydrogen atom.
I cannot tell the whole story here, but I can sketch the process.
When solving the Schrodinger equation involved, and transforming it to spherical coordinates, all applicable quantum numbers roll into your lap and on route you will encounter the Laguerre and associated Legendre differential equations.
Now for the approach of solving these equations. Let's have a look at the Laguerre equation first.
$$xy'' + (1+\alpha-x)y' + ny = 0$$
If you try $y = \sum_{j=0}^{\infty}a_jx^j$, you will find that
$$\frac{a_{j+1}}{a_j} = -\frac{n-j}{(j+1)(j+1+\alpha)}$$
So, when moving from one coefficient to the next, the sign changes, it is divided by $j+1$, multiplied by $n-j$ and divided by $j+1+\alpha$. When j reaches n, the coefficients vanish.
That means that we will not be far off if we try
$$a_j = \frac{(-1)^j}{j!}\binom{n+\alpha}{n-j}$$
In fact, this is exactly what we need.
Now, which function will this be? The factor $\frac{(-1)^j}{j!}$ points at $e^{-x}$.
The factor $\binom{n+\alpha}{n-j}$ points at the $j^{th}$ derivative of $x^{n+\alpha}$, and we saw above that the series ends for j=n. But when we try $[e^{-x}x^{n+\alpha}]^{(n)}$, firstly we have to compensate for all the $e^{-x}$ factors by multiplying the result by $e^x$ , and secondly we would get terms for the running index j that have a factor $x^{n + \alpha - (n-j)} = x^{j + \alpha}$ and we would like to have $x^j$. So we have to compensate at the end here too, by multiplying by $x^{-\alpha}$. This reasoning leads to trying
$$y = e^xx^{-\alpha}[e^{-x}x^{n+\alpha}]^{(n)}$$
Where $f^{(n)}$ stands for the $n^{th}$ derivative of $f$.
And funnily enough, it's spot on again.
Over to the Legendre equation, first the ordinary one.
$$[(1-x^2)y']' + n(n+1)y = 0$$
Trying $y = \sum_{j=0}^{\infty}a_jx^j$ again, we find that
$$\frac{a_{j+2}}{a_j} = \frac{j(j+1) - n(n+1)}{(j+1)(j+2)}$$
As can be seen in one of my other posts (About the Legendre differential equation) there are solutions that are polynomials ($P_n$) and solutions that are represented by an unending series expansion ($Q_n$).
The $P_n$ are either odd-powered or even powered polynomials, depending on $n$ which is also the degree of $P_n$.
The $Q_n$ contain a factor $ln(\frac{1+x}{1-x})$, see (About the Legendre differential equation) again.
Finally, the associated Legendre equation.
$$[(1-x^2)y']' + \{l(l+1) - \frac{m^2}{1-x^2}\}y = 0$$
It is natural to try $y = (1-x^2)^{\alpha}P_l$. It is easily checked that we get the desired term $\frac{m^2}{1-x^2}y$ if we take $\alpha = m/2$. Further analysis shows that if we wish to arrive at all the correct terms, we have to use $P_l^{(m)}$ instead of $P_l$. Note that I still use the notation $f^{(m)}$ as denoting the $m^{th}$ derivative of $f$. The notation I use for the solutions of the associated Legendre equation is
$$P_{l,m} = (1-x^2)^{m/2}P_l^{(m)}$$
I hope this helps :-).
Mathematically, not really.
I noticed that your formula for pendulum is not correct. You might mean these:
$$\frac{\dot{\theta}^2}{2} = \frac{g}{l} \cos(\theta)$$
and
$$\ddot{\theta} = -\frac{g}{l} \sin(\theta)$$
The trivial counterexample is $\theta \equiv \pi /2$, where $\dot{\theta}^2/2=\frac{g}{l}\cos\theta = 0$
but $\ddot{\theta}=0 \ne -\frac{g}{l}\sin\theta$.
In that case, you may assume that $\theta$ is not a constant.
By taking derivative, you get
$\ddot{\theta}\dot{\theta}=\frac{g}{l}\sin\theta \dot{\theta}$ and you can cancel $\dot{\theta}$ because $\theta$ is not a constant, then you get $\ddot{\theta} = -\frac{g}{l}\sin\theta$. However, if you have $\ddot{\theta} = -\frac{g}{l}\sin\theta$, then you can't get $\frac{\dot{\theta}^2}{2} = \frac{g}{l} \cos(\theta)$, because actually $\frac{\dot{\theta}^2}{2} = \frac{g}{l} \cos(\theta) + C$.
Best Answer
If the question is how do you get from your equation to the Legendre equation, from the context, I would expect you have $\theta \in [0,\pi]$. Then make the substitution $\cos \theta = x$, and the differential equation you have becomes, for $x \in[-1,1]$, \begin{align} (1-x^2) \frac {d^2 \Theta}{dx^2} - 2x \frac{d\Theta}{dx} + k\Theta = 0 \end{align} This is the Legendre equation. It has been studied extensively, and is described in the Wikipedia article referenced by Mr. Swanson in the comments. It has solutions known as the Legendre functions for general $k$, denoted $P_\lambda$ and $Q_\lambda$ where, by convention, $k=\lambda(\lambda+1)$ (I am using $\lambda$ instead of $\ell$ to match the article). The change of expression for $k$ is simply applied to turn your equation into a form that has recognized solutions.
When $\lambda$ is an integer, $P_\lambda$ a polynomial and the $Q_\lambda$ are known as a Legendre functions of the second kind. The boundary conditions for your equation will likely help determine which of these functions are relevant to your case. For example, applying boundary conditions may constrain $k$ to be an integer of the form $\lambda(\lambda+1)$ and you obtain the Legendre polynomials as an orthogonal family of solutions that forms a functional basis.