In here, I mentioned that the eigenvalues of a Hermitian operator are real and that the eigenfunctions of a Hermitian operator are orthogonal.

Let $\mathcal{L}$ be a Hermitian operator and let $u_i$, $u_j$ be eigenfunctions of $\mathcal{L}$ with eigenvalues $\lambda_i$, $\lambda_j$, respectively. Then
\begin{align}
\label{eq:eigen}
\mathcal{L}u_i+\lambda_iwu_i=0\\
\label{eq:eigen2}
\mathcal{L}u_j+\lambda_jwu_j=0
\end{align}

The complex conjugation of \eqref{eq:eigen2} is
\begin{equation}\label{eq:eigen3}\mathcal{L}u_j^\ast+\lambda_j^\ast wu_j^\ast=0\end{equation}
Multiply \eqref{eq:eigen} by $u_j^\ast$ and \eqref{eq:eigen3} by $u_i$:
\begin{align}
\label{eq:eigen4}u_j^\ast\mathcal{L}u_i+u_j^\ast\lambda_iwu_i=0\\
\label{eq:eigen5}u_i\mathcal{L}u_j^\ast+u_i\lambda_j^\ast wu_j^\ast=0
\end{align}
Subtracting \eqref{eq:eigen5} from \eqref{eq:eigen4}, we obtain
\begin{equation}\label{eq:eigen6}u_j^\ast\mathcal{L}u_i-u_i\mathcal{L}u_j^\ast=(\lambda_j^\ast-\lambda_i)wu_iu_j^\ast\end{equation}
Integrating \eqref{eq:eigen6}, we get
\begin{equation}\label{eq:eigen7}\int_a^bu_j^\ast\mathcal{L}u_i dx-\int_a^b u_i\mathcal{L}u_j^\ast dx=(\lambda_j^\ast-\lambda_i)\int_a^b u_iu_j^\ast wdx\end{equation}
Since $\mathcal{L}$ is Hermitian, the LHS of \eqref{eq:eigen7} vanishes. Thus
$$(\lambda_j^\ast-\lambda_i)\int_a^bu_iu_j^\ast wdx=0$$ If $i=j$ then $\int_a^bu_iu_j^\ast wdx\ne 0$ so $\lambda_i^\ast=\lambda_i$, i.e. $\lambda_i$ is real.

If $i\ne j$ and if $\lambda_i\ne\lambda_j$, then $\int_a^bu_iu_j^\ast wdx=0$, i.e. $u_i$ and $u_j$ are orthogonal, where we consider
\begin{equation}\label{eq:hermitianproduct}\langle u_i,u_j\rangle=\int_a^b u_iu_j^\ast wdx\end{equation} as an inner product. The inner product \eqref{eq:hermitianproduct} is called the Hermitian product weight by $w(x)$.

What if $i\ne j$ but $\lambda_i=\lambda_j$? Then $\int_a^bu_iu_j^\ast wdx$ may not vanish, i.e. $u_i$ and $u_j$ may not be orthogonal. Such case is labeled degenerate. For example, consider the differential equation
$$\frac{d^2}{dx^2}y(x)+n^2y(x)=0$$
Once can easily see that $y_1(x)=\cos nx$ and $y_2(x)=\sin nx$ both satisfy the differential equation. So $\cos nx$ and $\sin nx$ are eigenfunctions that correspond to the same eigenvalue $n^2$. That is, $\cos nx$ and $\sin nx$ are degenerate eigenfunctions. They are however orthogonal because
$$\int_0^{2\pi}\cos nx\sin nxdx=0$$

It is good to know the following formulas:
$$
\int_{x_0}^{x_0+2\pi}\sin mx\sin nxdx=C_n\delta_{nm},
$$
where
$$C_n=\left\{\begin{array}{ccc}
\pi & \mbox{if} & n\ne 0\\
0 & \mbox{if} & n=0
\end{array}\right.$$
$$
\int_{x_0}^{x_0+2\pi}\cos mx\cos nxdx=D_n\delta_{nm},
$$
where
$$D_n=\left\{\begin{array}{ccc}
\pi & \mbox{if} & n\ne 0\\
2\pi & \mbox{if} & n=0
\end{array}\right.$$
$$\int_{x_0}^{x_0+2\pi}\sin mx\cos nxdx=0$$

Orthogonal functions can be used to expand a functions, for example, as a Fourier series expansion. Certain classes of function (sectionally continuous or piecewise continuous) may be represented by a series of orthogonal functions to any desired degree of accuracy. Such property is referred to as completeness.

Example. Square Wave

Consider the square wave
$$f(x)=\left\{\begin{array}{ccc}
\frac{h}{2} & \mbox{if} & 0<x<\pi,\\
-\frac{h}{2} & \mbox{if} & -\pi<x<0
\end{array}\right.$$

Note that the function may be expanded in any of a variety of orthogonal eigenfunctions such as Legendere, Hermite, Chebyshev, etc. The choice of eigenfunction is made on the basis of convenience. For example, we may choose to use $\cos nx$, $\sin nx$.  The eigenfunction series can then be written as a Fourier series
\begin{equation}\label{eq:fourier}f(x)=\frac{a_0}{2}+\sum_{n=1}^\infty(a_n\cos nx+b_n\sin nx)\end{equation}
where
\begin{align*}
a_n&=\frac{1}{\pi}\int_{-\pi}^\pi f(x)\cos nt dt,\\
b_n&=\frac{1}{\pi}\int_{-\pi}^\pi f(t)\sin nt dt,\ n=0,1,2,\cdots
\end{align*}
So by the formula \eqref{eq:fourier} we find that
\begin{align*}
a_n&=0,\\
b_n&=\frac{h}{n\pi}(1-\cos n\pi)\\
&=\frac{h}{n\pi}(1-(-1)^n)\\
&=\left\{\begin{array}{ccc}
0 & \mbox{if} & n=\mbox{even}\\
\frac{2h}{n\pi} & \mbox{if} & n=\mbox{odd}
\end{array}\right.
\end{align*}
Hence the Fourier series expansion of $f(x)$ is given by
$$f(x)=\frac{2h}{\pi}\sum_{n=0}^\infty\frac{\sin(2n+1)x}{2n+1}$$

If $N$ liearly independent eigenfunctions correspond to the same eigenvalue, the eigenvalue is said to be $N$-fold degenerate. In the above example, for each $n$ there are two possible solutions $\cos nx$, $\sin nx$. We may say the eigenfunctions are degenerate or the eigenvalue is degenerate.

Gram-Schmidt Orthogonalization

Consider
$$\int_a^b \varphi_i^2wdx=N_i^2$$
Since $\mathcal{L}u(x)+\lambda w(x)u(x)=0$ is linear, $\mu_i\varphi_i$ would be a solution as well for any constant $\mu_i$. If we set $\psi_i=\frac{\varphi_i}{N_i}$ then
$$\int_a^b\psi_i^2wdx=1$$
and
$$\int_a^b\psi_i(x)\psi_j(x)w(x)dx=\delta_{ij}$$

Suppose that $\{u_n(x)\}_{n=0}^\infty$ are eigenfunctions that are not mutually orthogonal. Gram-Schmidt orthogonalization process allows us to come up with orthonormal eigenfunctions $\{\varphi_n\}_{n=0}^\infty$ from $\{u_n(x)\}_{n=0}^\infty$.
Let us start with $n=0$. Let $\psi_0(x)=u_0(x)$. Then the normalized eigenfunction is denoted by $\varphi_0(x)$.
$$\varphi_0(x)=\frac{\psi_0(x)}{\left[\int\psi_0^2(x)wdx\right]^{1/2}}$$
For $n=1$, let
\begin{align*}
\psi_1&=u_1(x)-\langle u_1(x),\varphi_0(x)\rangle\varphi_0(x)\\
&=u_1(x)-\int u_1\varphi_0wdx\varphi_0(x)\end{align*}
and
$$\varphi_1(x)=\frac{\psi_1}{\left[\int\psi_1^2wdx\right]^{1/2}}$$

Figure 1

Figure 1 clearly shows how we can come up with $\psi_1$ using a vector projection. Note that the angle $\theta$ appeared in Figure 1 is for an intuitive purpose only and does not depict a real angle.

Figure 2

For $n=2$, as one can see from Figure 2, $\psi_2$ that is orthogonal to both $\varphi_0$ and $\varphi_1$ can be obtained as
\begin{align*}
\psi_2&:=u_2(x)-\langle u_2(x),\varphi_0(x)\rangle\varphi_0(x)-\langle u_2(x),\varphi_1(x)\rangle\varphi_1(x)\\
&=u_2(x)-\int u_2(x)\varphi_0(x)w(x)dx\varphi_0(x)-\int u_2(x)\varphi_1(x)w(x)dx\varphi_1(x)
\end{align*}
The normalzed eigenfunction $\varphi_2$ is then
$$\varphi_2:=\frac{\psi_2}{\left[\int\psi_2^2wdx\right]^{1/2}}$$
Now we can see a clear pattern and $\psi_n$ which is orthogonal to $\varphi_0,\cdots,\varphi_{n-1}$ would be given by
\begin{align*}
\psi_n&:=u_n(x)-\sum_{j=0}^{n-1}\langle u_j(x),\varphi_j(x)\rangle\varphi_j(x)\\
&=u_n(x)-\sum_{j=0}^{n-1}\left[\int u_j(x)\varphi_j(x)w(x)dx\right]\varphi_j(x)
\end{align*}
with its normalization
$$\varphi_n:=\frac{\varphi_n}{\left[\int\psi_n^2wdx\right]^{1/2}}$$

Example. Find orthonormal set from
$$u_n(x)=x^n,\ n=0,1,2,\cdots$$
with $-1\leq x\leq 1$ and $w(x)=1$.

Solution. By Gram-Schmidt orthonormalization, we obtain
\begin{align*}
\varphi_0(x)&=\frac{1}{\sqrt{x}},\\
\varphi_1(x)&=\sqrt{\frac{3}{2}}x,\\
\varphi_2(x)&=\sqrt{\frac{5}{2}}\cdot\frac{1}{2}(3x^2-1),\\ \varphi_3(x)&=\sqrt{\frac{7}{2}}\cdot\frac{1}{2}(5x^3-3x),\\
\vdots\\
\varphi_n(x)&=\sqrt{\frac{2n+1}{2}}P_n(x),
\end{align*}
where $P_n(x)$ is the $n$th-order Legendre polynomial.

References:

G. Arfken, Mathematical Methods for Physicists, 3rd Edition, Academic Press 1985