Since $\phi_n$ and $\phi_m$ are eigenfunctions, they must satisfy the ODE ... I came across this:

Theorem 1
The eigenfunctions of Sturm-Liouville BVP above satisfy the integral relationship:
$$\int_a^b r(x)\phi_n(x) \phi_m(x) \ dx = 0$$ if $m \not= n$,
where $\phi_1, \phi_2, \phi_3, \dots$ are eigenfunctions and $\phi_n$ corresponds to the eigenvalue $\lambda_n$.
Hence, the set of eigenfunctions for the Sturm-Liouville problem are orthogonal on the interval of interest w.r.t the weight function $r(x)$.

Then, when proving this, it starts with

Since $\phi_n$ and $\phi_m$ are eigenfunctions, they must satisfy the ODE
$$\frac{d}{dx}\left(p(x) \frac{d \phi_n}{dx}\right) + q(x) \phi_n = - \lambda_n r(x) \phi_n$$
$$\frac{d}{dx}\left(p(x) \frac{d \phi_m}{dx}\right) + q(x) \phi_m = - \lambda_m r(x) \phi_m$$

I was stumped by this. I'm wondering why, since $\phi_n$ and $\phi_m$ are eigenfunctions, they must satisfy those ODEs?
 A: Sturm-Liouville equations originally arose out Fourier's separation of variables method of solving partial differential equations. The separation parameter $\lambda$ serves the role of an eigenvalue for an eigenfunction problem. In fact, eigenvector/eigenvalue analysis of matrices came out of studying Sturm-Liouville equations using Fourier's method. Sturm and Liouville isolate a general type of equation that would arise out of Fourier separation problems, and they set out of study this equation, and to show that every function $f$ on the interval of interest $[a,b]$ could be expanded in a Fourier series of eigenfunctions of the Sturm-Liouville problem:
$$
             -\frac{d}{dx}\left(p(x)\frac{df}{dx}\right)+q(x)f(x)= \lambda w(x)f(x),\\
        A_1 f(a)+A_2 f'(a) = 0,\;\;\; B_1 f(b)+B_2 f'(b)=0.
$$
For a regular problem where $p$ does not vanish at $a$ or $b$, there are only a discrete number of values $\lambda_1,\lambda_2,\lambda_3,\cdots$ for which there are corresponding solutions $f_n$ that are not identically $0$ and satisfy both endpoint conditions. These are the eigenvalues, and $f_n$ is the eigenfunction that corresponds to $\lambda_n$. For these solutions, it is possible to write any smooth function $g$ on $[a,b]$ in terms of these functions in a Fourier series $g=\sum_{n=1}^{\infty}c_n f_n(x)$ where $c_n$ are constants. The constants are determined by the orthogonality conditions that arise from the equation:
$$
       \int_{a}^{b}f_n(x)f_m(x)w(x)dx = 0,\;\;\; n\ne m.
$$
In fact, assuming $g=\sum_{n=1}^{\infty}c_n f_n(x)$, you multiply by $f_m(x)w(x)$, integrate, and use the above to obtain the unknown constants $c_m$:
$$
         \int_a^b g(x)f_m(x)w(x)dx = \sum_{n=1}^{\infty}c_n \int_{a}^{b}f_m(x)f_n(x)w(x)dx \\
       \int_a^b g(x)f_m(x)w(x)dx = c_m\int_a^b f_m(x)^2w(x)dx \\
             c_m  = \frac{\int_a^b g(x)f_m(x) w(x)dx}{\int_a^b f_m(x)^2w(x)dx}
$$
And, in a fairly general sense, you can expand a function $f$ in a Fourier series of the eigenfunctions that has the form
$$
         g(x) = \sum_{m=1}^{\infty}\frac{\int_a^b g(y)f_m(y) w(y)dy}{\int_a^b f(y)_m^2w(y)dy} f_m(x).
$$
Finding the eignevalues, the eigenfunctions, and expanding in such a Fourier series is the goal of Sturm-Liouville Theory. Using this, you can solve various PDEs using Fourier's separation of variables technique.
