A while back I found a specific proof for the orthogonality of two functions but I can't seem to find it online anywhere. I just want to make sure the proof I've given is defined exactly, which is:
"Now more formally take the two function $\phi_m, \phi_n$ and they are orthogonal under some interval $[a,b]$. Therefore, \begin{align*} \langle \phi_n, \phi_m \rangle & = \int_a^b \phi_m \cdot \phi_n dx \\ & \to 0 \Longleftrightarrow n \neq m \\ & \to ||\phi_n^2|| \Longleftrightarrow n = m \end{align*} Where $|| \vec{v} ||$ is the square norm of a vector, which is just an absolute value but for vectors so that they're positive in all directions. Now take another function $f(x)$, such that $$f(x) = \sum_{n=0}^{\infty} C_n \phi_n $$ $$ \langle f(x), \phi_m \rangle = \sum_{n=0}^{\infty} C_n\langle \phi_n, \phi_m \rangle $$
$$ \langle f(x), \phi_m \rangle = \sum_{n=0}^{\infty} C_n \int_a^b \phi_m \cdot \phi_n dx $$ Which can be simplified using the aforementioned property, $$ \langle \phi_n, \phi_m \rangle = 0 \Longleftrightarrow n \neq m $$ Now $C_n$, the coefficient, can be computed.
\begin{align*} \langle f, \phi_m \rangle & = C_n\langle \phi_m, \phi_m \rangle \\ C_n & = \frac{\langle f, \phi_m \rangle}{\langle \phi_m, \phi_m \rangle} \\ & = \frac{\langle f, \phi_m \rangle}{||\phi_n^2||} \end{align*} This principle is then used with sinusoidal functions throughout this chapter, and it is this key property that underlies all of Fourier Series and Transforms."
It comes about from a Sturm Liouville eigenvalue problem
$$ \frac{d}{dx}\bigg( p(x) \frac{d\phi}{dx}\bigg) + q(x)\phi + \lambda\sigma(x) \phi =0 \tag{1} $$
where $\lambda$ is an eigenvalue and $x$ is on an interval $a < x < b$ Eigenfunctions belonging to different eigenvalues are orthognal with respect to a weight function $\sigma(x)$
$$ \int_{a}^{b} \phi_{n}(x)\phi_{m}(x)\sigma(x) dx = 0 , \lambda_{n} \neq \lambda_{m} \tag{2}$$
the most basic case we can show is for what you said. Let $p(x) = 1 q(x) =0 \sigma(x) =1,$ then the first equation becomes
$$ \frac{d^{2}\phi}{dx^{2}} +\lambda\phi =0 \tag{3} $$
if we let the boundaries be
$$ \phi(0) = 0 \\ \phi(L) = 0 \tag{4} $$
and you get
$$ \lambda = \bigg(\frac{n\pi}{L}\bigg)^{2} \tag{5}$$
$$ \phi(x) = \sin(\sqrt{\lambda}x) = \sin(\frac{n\pi x}{L}) \tag{6} $$
The weight function as we noted is $\sigma(x) = 1$
$$ \int_{0}^{L} \sin(\frac{n\pi x}{L})\sin(\frac{m\pi x}{L}) dx =0 \tag{7}$$
In general, you use something called Fourier's trick. I don't think really answers your question but it doesn note that eigenfunctions aren't orthogonal in general. They need a weight function $\sigma(x)$