My teacher talked about Sturm-Liouville theory, and we learned that any second order differential equation can be put into the self-adjoint form. What is this? Well, the book says is something in the form:
$$\frac{d}{dx}\left[p(x)\frac{du(x)}{dx}\right] + q(x)u(x) \tag {1}$$
The book also says that we can put any second order linear ODE into the self-adjoint form:
Consider the second order linear ODE:
$$Lu(x) = p_0\frac{d^2}{dx^2}u(x) + p_1(x)\frac{d}{dx}u(x) + p_2(x)u(x)$$
Multiply it by:
$$\frac{1}{p_0(x)}\exp\left[\int^x\frac{p_1(t)}{p_0(t)}\ dt\right]$$
to obtain:
$$\frac{1}{p_0(x)}\exp\left[\int^x\frac{p_1(t)}{p_0(t)}\ dt\right]Lu(x) = \frac{d}{dx}\left\{\exp\left[\int^x\frac{p_1(t)}{p_0(t)}\ dt\right]\frac{du(x)}{dx}\right\} + \frac{p_2(x)}{p_0(x)}\cdot \exp\left[\int^x\frac{p_1(t)}{p_0(t)} \ dt\right]u$$
Which is in the self-adjoint form, whatever that means. We can see by comparing with eq $(1)$ and taking $p(x) = \exp\left[\int^x\frac{p_1(t)}{p_0(t)}\ dt\right]$ and $q(x) = \frac{p_2(x)}{p_0(x)}\cdot \exp\left[\int^x\frac{p_1(t)}{p_0(t)} \ dt\right]$.
Then the book, and also my teacher, says that we're interested in finding eigenfunctions and eigenvalues for:
$$\frac{d}{dx}\left[p(x)\frac{du(x)}{dx}\right] + q(x)u(x) + \lambda w(x) u(x) = \\ Lu(x) + \lambda w(x) u(x) $$
for some density function $w(x)$ (what), where $\lambda$ is a eigenvalue related to the eigenfunction $u(x)$, called in the book $u_{\lambda}$ because it depends on $\lambda$.
The book and my teacher then talk about completeness of eigenvalues and eigenfunctions, and that's it. Nothing more is said about Sturm-Liouville theory, then we jumped to other things and I'm asking for what all this is used.
It also mentions boundary conditions. What are boundary conditions? Are them simply the initial conditions for some ODE? Because the boundary conditions given in the book are of the form
$$a_1 u(a) + a_2 u'(a) = 0\\ b_1 u(a) + b_2 u'(a) = 0$$
but aren't initial conditions things in the form $u'(a) = b, u''(a) = c$?
When I see examples on the wikipedia, they all mention that we already know that such eigenfunction $u(x)$ is a solution with eigenvalue $\lambda$, and also in these slides I see no examples of how to actually finding the eigenvalue and eigenfunction...
So what is the usefulness of the self-adjoint formula and of the sturm-liouville theory in general?
Sturm-Liouville Theory came out of studying Fourier's Separation of Variables technique for solving Partial Differential Equations, and that's still an important application. In this context, solving the ODE is not the ultimate objective, but it provides a way to solve a more complicated equation. As an example, suppose $A$ is a selfadjoint matrix, and you want to solve a vector equation $$ \frac{d x}{dt} = A x,\;\;\; x(0)=x_0. $$ This problem can be reduced by finding an orthonormal basis $\{ e_1,e_2,\cdots,e_n \}$ of eigenvectors for $A$ with eigenvalues $\lambda_1,\lambda_2,\cdots\lambda_n$ because you can then write $$ x(t) = \alpha_1(t)e_1 + \cdots + \alpha_n(t)e_n $$ where $\alpha_k(t)$ are scalar functions of $t$, and plug back into the equation to obtain $$ \sum_{k=1}^{n} \alpha_k'(t)e_k = \sum_{k=1}^{n} \lambda_k \alpha_k(t)e_k \\ \implies \alpha_k(t) = \alpha_k(0)e^{\lambda_k t}. $$ Then you can find $\alpha_k(0)$ knowing that $x(0)=x_0$: $$ x_0 = x(0) = \sum_{k=0}^{n}\alpha_k(0)e_k \\ \implies \langle x_0,e_l\rangle = \alpha_l(0) \\ x(t) = \sum_{k=1}^{n}e^{\lambda_k t}\langle x_0,e_k\rangle e_k $$ This is a full solution of the problem.
Separation of variables works basically the same way for partial differential equations, but on a grander scale. Here Sturm-Liouville operators must be diagonalized, which requires an infinite-dimensional space, and an infinite orthogonal basis. The operators are selafadjoint, and so there is a strong similarity between these cases and the matrix case mentioned above. For example, the wave equation for a thin uniform wire stretched between two fixed points and plucked to an initial displacement $u_0(x)$ initially at rest and then released at $t=0$ is a partial differential equation for the unknown displacement $u(t,x)$ of the string governed by $$ \frac{\partial^2 u}{\partial t^2}=c\frac{\partial^2 u}{\partial x^2},\\ u(t,0)=u(t,L)=0, \\ u(0,x)=f(x), \\ u_{t}(0,x)=0. $$ In this case you're working with an operator $A f = -\frac{\partial^2}{\partial x^2}f$ on the domain of twice differentiable functions $f$ for which $f(0)=f(L)=0$. This operator also has an orthonormal basis of eigenfunctions, just as a selfadjoint matrix does: $$ e_n(x) = \sqrt{\frac{2}{L}}\sin(n\pi x/L). $$ The operator has eigenvalues $\lambda_n = n^2\pi^2/L^2$ and $$ Le_n = \lambda_n e_n. $$ By expanding $u(t,x)$ in these functions, the $x$ variable is removed, and you end up with a time-varying vector equation where the vectors are functions of $x$ in the vector space $L^2[0,L]$. $$ u(t) = \sum_{n=1}^{\infty}\alpha_n(t)e_n $$ The coefficients $\alpha_n(t)$ are scalar coefficient functions that must satisfy $$ \alpha_n''(t) = c \alpha_n(t),\;\;\; \alpha_n(0) = \langle f,e_n\rangle,\;\; \alpha_n'(0)=0. $$ The coefficient function solutions are $$ \alpha_n(t) = \langle f,e_n\rangle\cos(\sqrt{c}n\pi t/L), $$ which leads to the PDE solution $$ u(t,x) = \sum_{n=1}^{\infty}\langle f,e_n\rangle \cos(\sqrt{c}n\pi t/L)e_n(x) \\ = \sum_{n=1}^{\infty}\left[\frac{2}{L}\int_{0}^{L}f(x')\sin(n\pi x'/L)dx'\right]\cos(\sqrt{c}n\pi t/L)\sin(n\pi x/L) $$ This solution is a decomposition of the vibrational discplacement in terms of the "natural modes" of the string, which are harmonic multiples of a fixed ground state frequency, which is how the field came to be known as Harmonic Analysis.
So the goal is not so much how to solve the Sturm-Liouville problem, but how to use the "diagonalization" of the Sturm-Liouville problem in terms of its eigenfunctions to reduce partial differential equation problems to ODE problems for the coefficients in eigenfunction expansions of the solution along a basis that diagonalizes a Sturm-Liouville operator.