Suppose you are working on a bounded interval $[a,b]$ with Sturm-Liouville operator $L$ given by
$$
Lf = \frac{1}{w(x)}\left[-\frac{d}{dx}\left(p(x)\frac{df}{dx}\right)+q(x)f\right].
$$
Assume that $p,q,w$ are real functions, that $p$ is continuously differentiable, and that $q$, $w$ are continuous on $[a,b]$ with $w > 0$ and $p > 0$ on $[a,b]$. The Sturm-Liouville eigenvalue problem for $L$ is posed with the domain $\mathcal{D}(L)$ restricted to twice continuously differentiable funtions $f$ on $[a,b]$ that satisfy endpoint conditions
$$
\cos\alpha f(a)+\sin\alpha f'(a) =0,\\
\cos\beta f(b) + \sin\beta f'(b) = 0.
$$
The natural space to consider this problem is in $L^2_{w}[a,b]$, with inner product
$$
(f,g)_w = \int_{a}^{b}f(x)\overline{g(x)}w(x)dx.
$$
For $f,g\in\mathcal{D}(L)$ (meaning both satisfy the endpoint conditions,) the Lagrange identity gives
$$
(Lf,g)_w = (f,Lg)_w.
$$
Any time you have a symmetric operator on a complex inner product space, the eigenvalues must be real. To see why, suppose $Lf = \lambda f$. Then
$$
\lambda (f,f)_w = (Lf,f)_w = (f,Lf)_w = (f,\lambda f)_w = \overline{\lambda}(f,f)_w.
$$
So, either (a) $\lambda=\overline{\lambda}$ or (b) $(f,f)_w=0$, which means $f=0$. The operator $L$ is real, meaning that $\overline{Lf}=L\overline{f}$. Therefore, if $Lf=\lambda f$, then $\lambda$ must be real, and
$$
L\overline{f} = \overline{Lf} = \overline{\lambda f}=\lambda\overline{f}.
$$
In other words, if $f$ is an eigenfunction, then so is $\overline{f}$, which means that $\Re f$ and $\Im f$ are eigenfunctions, and both $f$, $\overline{f}$ satisfy the same endpoint conditions because the endpoint conditions are formulated with real coefficients. To see that $\{ \Re f,\Im f\}$ must be a linearly dependent set of solutions, start by noticing that if $g$ is a real eigenfunction with $g(a)=g'(a)=0$, then $g$ must be zero by classical uniqueness results for second order ODE's. Therefore, $\Re f$ must satisfy the following for some non-zero scalar $C$:
$$
(\Re f)(a) = -C\sin\alpha,\;\; (\Re f)'(a)=C\cos\alpha
$$
And the same must be true of $\Im f$, with a different constant $D \ne 0$. It follows that $g=D\Re f - C\Im f$ is a real eigenfunction with $g(a)=g'(a)=0$, which forces $g\equiv 0$. So $\{ \Re f,\Im f\}$ must be a linearly dependent set of functions. The final conclusion is that $f= \rho\Re f$ where $\rho$ is a non-zero complex number, or $f=\rho\Im f$ where $\rho$ is a non-zero complex number. So you can always assume the eigenfunctions are real after multiplying by the correct scalar.
Best Answer
It may be helpful to understand the similar fact about left and right eigenvectors of a matrix, in order to appreciate the role being played here by self-adjointness of the operator $L$.
In any case the computation is straightforward, applying what we know in the case that $u,v$ are eigenfunctions (satisfying the unspecified but requisite homogeneous boundary conditions) corresponding to distinct eigenvalues $\lambda_1 \neq \lambda_2$ respectively:
$$ L(u) = -\lambda_1 \sigma u $$
$$ L(v) = -\lambda_2 \sigma v $$
where multiplier function $\sigma(x)$ on $[a,b]$ must be assumed positive (a.e.) there if we want the natural bilinear function to be an inner product, with the corresponding sense of "orthogonality":
$$ 0 = \int_a^b u L(v) - v L(u) \mathrm{d} x = (\lambda_1 - \lambda_2) \int_a^b u v \; \sigma(x) \mathrm{d} x$$
Since the eigenvalues are unequal, the last integral is zero, and thus $u,v$ are orthogonal with respect to this weighted inner product on $[a,b]$.