Here are two reasons why Sturm-Liouville theory is useful:
Physics reason: I noticed that you tagged this question as "physics". Almost every problem in quantum physics is a Sturm-Liouville eigenfunction problem: namely, solving Schrodinger's equation!
For example:
Take the infinite square well. Here, Schrodinger's equation is of the form
$ \mathcal L u + \lambda u = 0$, where $\mathcal L = \frac {d^2}{dx^2}$ (up to a constant factor). The boundary conditions $u(0) = u(\pi) = 0$. This is a Sturm-Liouville eigenfunction problem. The eigenfunctions $u_n = \sin(n x)$ are the wavefunctions of the various states, and the eigenvalues $\lambda_n = n^2$ are their energies. (The allowed values of $n$ are $n = 1,2,3,\dots$)
For the simple harmonic oscillator, Schrodinger's equation can be converted to the Hermite equation, which is also of Sturm-Liouville form. The eigenfunctions and eigenvalues are (closely related to) the wavefunctions and energies of the simple harmonic oscillator states. (This example is a "singular" Sturm-Liouville system, because the domain extends to infinity.)
When you solve Schrodinger's equation for the hydrogen atom by separation of variables, the radial equation can be converted into the associated Laguerre's equation, which is another (singular) Sturm-Liouville system...
Sturm-Liouville theory tells us (with certain caveats relating to dimensionality, infinite domains and singularities, and with subtleties depending on how closely-related the Sturm-Liouville problem is to the original Schrodinger problem) that:
There is a discrete sequence of energy levels, labelled by a quantum number $n$. This discreteness is why quantum physics is called "quantum". The energy values of the states form an ascending sequence $\lambda_1 \leq \lambda_2 \leq \lambda_3 \leq \dots$, which tends to infinity. So in particular, there is a well-defined notation of a "ground state" of minimal energy: this is the $n = 1$ state.
The wavefunctions of two energy levels with distinct energy values are orthogonal. This is an easy consequence of the self-adjointness property of the Hamiltonian operator. If we rescale the wavefunctions (and choose bases carefully if certain energy levels are degenerate), then we can make the wavefunctions orthonormal: i.e. $\int u_{n_1} u_{n_2} = \delta_{n_1, n_2}$.
Any general wavefunction $u$ can be written as a linear combination $ u = \sum_n c_n u_{n} $ of the pure energy-level wavefunctions $u_{n}$. This is the "completeness" statement in Sturm-Liouville theory. In physics language, the state $u$ is then a "superposition" of the pure energy-level states.
Moreover, it is easy to work out the coefficients $c_n$ in the linear decomposition $ u = \sum_n c_n u_{n} $. Since we have an orthogonality relation $\int u_{n_1} u_{n_2} = \delta_{n_1 , n_2}$, the $c_n$ coefficient is simply $c_n = \int u_{n} u$.
Maths reason: Suppose we wish to solve an equation of the form
$$ \mathcal L u = f,$$
where $\mathcal L$ is a second-order differential operator of Sturm-Liouville type, and $f$ is some given function. (For example, this could be a dynamical system, where $f$ represents a forcing term.)
A common strategy for solving an equation of this form is to start by finding a set of eigenvalues $\lambda_n$ and eigenfunctions $u_n$ satisfying the equation,
$$ \mathcal L u_n + \lambda_n w u_n = 0,$$
You're right in saying that Sturm-Liouville theory doesn't actually help you to find these eigenvalues and eigenfunctions. Usually, people determine them using other methods, such as power series. [Many of these power series solutions to Sturm-Liouville eigenfunction problems are well known - perhaps you can read up on Legendre polynomials, Hermite polynomials, Laguerre polynomials, Chebyshev polynomials, and so on.]
But once we have found the $\lambda_n$'s and $u_n$ satisfying the eigenfunction equation, we can use these eigenvalues and eigenfunctions, together with results from Sturm-Liouville theory, to build solutions to the original equation $\mathcal L u = f$. Here is how we do it:
First, Sturm-Liouville theory tells us that the $u_n$'s are complete. So we can write the solution $u$ of our original equation as a linear combination of the $u_n$'s:
$$ u = \sum_n c_n u_n.$$
If we now substitute this ansatz into the equation $\mathcal L u = f$, we get
$$- \sum_n c_n \lambda_n w u_n = f$$
Sturm-Liouville theory also tells us that the $u_n$'s are orthonormal (after rescaling): $$\int w u_{n_1} u_{n_2} = \delta_{n_1, n_2}.$$
[By a standard integration-by-parts argument, this orthonormality property follows from the fact that a differential operator of the form $\mathcal L = \frac d {dx} p \frac d {dx} + q$ is self-adjoint.]
So if we multiply both sides by $u_k $ and integrate, we get
$$ c_k = - \tfrac 1 {\lambda_k}\int f u_k .$$
We have now solved the equation! The answer is
$$ u = \sum_n \left( - \tfrac 1 {\lambda_n}\int f u_n \right) u_n.$$
[If you like, you can think of $\sum_n \tfrac 1 {\lambda_n }u_n(x')u_n(x)$ as the Green's function for the operator $\mathcal L$.]
To summarise, Sturm-Liouville theory doesn't tell us how to find the eigenfunctions and eigenvalues that satisfy $\mathcal L u + \lambda w u = 0$; this is done usually using power series methods. Instead, Sturm-Liouville theory tells us information about these eigenfunctions and eigenvalues that enable us to use them as building blocks for building solutions to general problems of the form $\mathcal L u = f$.
[This technique can be generalised for equations of the form
$ \mathcal L u + \mu w u = f,$
where $\mu$ is a constant number. As along as $\mu$ is not equal to any $\lambda_n$, one can use the same method that the solution is
$ u = \sum_n \left( \tfrac 1 {\mu - \lambda_n}\int f u_n \right) u_n.$
Since Sturm-Liouville theory tells us that the $\lambda_n$'s form a discrete set (more specifically, a positive, strictly-increasing sequence), it follows that this solution method is valid for "almost" every choice of $\mu$.]
There are other substitutions you can use to transform the equation to a normal form. For example, the equation
$$
a(x)y''(x)+b(x)y'(x)+c(x)y(x)+\lambda d(x)y(x)=0,
$$
is transformed to
$$
-\frac{d}{dx}\left(a(x)\frac{df}{dx}\right)+\left(\frac{(b-a')^2}{4a}+\frac{(b-a')'}{2}-c\right)f = \lambda d(x)f(x)
$$
by the substitution $y=\rho(x)f(x)$, where
$$
\rho(x) = \sqrt{a(x)}\exp\left(-\int\frac{b}{2a}dx\right).
$$
This substitution is far more versatile in dealing with classical equations. However, it requires $a(x)$ to be twice differentiable. The advantage is that $a(x)$ remains the coefficient of $f''$ and $d(x)$ remains the of the eigenvalue $\lambda f(x)$, which is what I believe you're after. Many classical standard forms can be derived using this substitution.
For example, the associated Legendre equation
$$
(1-x^2)u''-2x(m+1)u'+[l(l+1)-m^2-m]u=0
$$
is transformed by $u=(1-x^2)^{-m/2}f$ to the classical form
$$
-\frac{d}{dx}\left((1-x^2)\frac{df}{d}\right)+\frac{m^2}{1-x^2}f=l(l+1)f.
$$
This form is hard to get any other way.
Another use of this substitution is to transform an equation to potential form
$$
-\frac{d^2f}{dx^2}+qf = \lambda f.
$$
This requires a clever use of the above substitution, followed by a substitution of the independent variable.
ANOTHER APPROACH: Typically, the function $d$ is turned into a weight function for the inner product space when you have a form such as the one I mentioned above:
$$
-\frac{d}{dx}\left(a(x)\frac{df}{dx}\right)+\left(\frac{(b-a')^2}{4a}+\frac{(b-a')'}{2}-c\right)f = \lambda d(x)f(x)
$$
Then you can define the operator
$$
L=\frac{1}{d}\left[ -\frac{d}{dx}\left(a(x)\frac{df}{dx}\right)+\left(\frac{(b-a')^2}{4a}+\frac{(b-a')'}{2}-c\right)f \right]
$$
on the weighted $L^2$ space $L^2_{d}$ defined by
$$
\langle f,g\rangle_{L^2_{d}}=\int f(x)\overline{g(x)}d(x)dx
$$
In this setting the operator $L$ is symmetric with the right endpoint conditions:
$$
\langle Lf,g\rangle_{L^2_d}=\langle f,Lg\rangle_{L^2_d}+\mbox{evaluation terms},
$$
and the eigenvalue problem is $Lf=\lambda f$.
Best Answer
Following up on my comment, given $$y'' + 2y' + (\lambda + 1)y = 0$$ You can multiply by the "integrating factor" $\mu(x) = e^{\int 2 dx} = e^{2x}$, $$e^{2x}y'' + 2e^{2x}y' + (\lambda + 1)e^{2x}y = 0$$ At which point you have $$ \frac{d}{dx}\left(e^{2x} y' \right) + (e^{2x} + \lambda e^{2x})y = 0 $$ as desired.