A Hamiltonian is invariant under time reversals if $[H , T ] = 0$. Remember that $T$ is an anti-unitary operator. Now, consider the state $|n\rangle$ which satisfies
$$
H|n\rangle = E_n |n\rangle
$$
Now, consider the state $T|n\rangle$. Since $H$ commutes with $T$ and $E_n$ is real, it can easily be shown that the state $T|n\rangle$ also has energy $E_n$ (check this yourself). Now, if the state $|n\rangle$ is non-degenerate, then we must have
$$
T |n\rangle = e^{i\theta_n}|n\rangle , \qquad \theta_n \in [0,2\pi) \tag{1}
$$
$\theta_n$ is some phase. Now, lets look at the wave-function. Act on both sides of the equation above with $\langle x|$. Since time reversals does not act on the operator $x$, we have $T| x \rangle = | x \rangle$. Thus,
$$
\langle x | T |n\rangle = e^{i\theta_n}\langle x |n\rangle \implies e^{i\theta_n} \psi_n(x) = \langle n | T^\dagger|x\rangle = \langle n | T^{-1} |x\rangle = \langle n |x\rangle = \langle x |n \rangle ^* = \psi_n(x)^*
$$
Here, we have used the fact that for anti-linear operators, the adjoint is defined as $\langle n | T^\dagger x\rangle = \langle x | T n\rangle$ (instead of $\langle n | L^\dagger x\rangle = \langle L n | x \rangle$ which is the definition for linear operators). We have also used that the time-reversal operator is unitary so $T^\dagger = T^{-1}$.
The final result is
$$
\psi_n(x)^* = e^{i\theta_n} \psi_n(x). \tag{2}
$$
Here, $\theta_n$ does not itself depend on $x$. We now write
$$
\psi_n(x) = e^{-i\theta_n/2} f_n(x)
$$
Substituting this into equation (2), we find
$$
[ e^{-i\theta_n/2} f_n(x) ] ^* = e^{i\theta_n} [ e^{-i\theta_n/2} f_n(x) ] \implies e^{i\theta_n/2} f_n(x)^* = e^{i\theta_n} e^{-i\theta_n/2} f_n(x) \implies f_n(x)^* = f_n(x) .
$$
Thus, we find that we can write
$$
\psi_n(x) = e^{-i\theta_n/2} f_n(x) \quad \text{where} \quad f_n(x)^* = f_n(x) .
$$
QED.
Let me also discuss what happens if we relax the non-degeneracy requirement. In this case, there is no relation between the states $|n\rangle$ and $T|n\rangle$. Now recall that for all states in the Hilbert space, $T^2$ acts as $+1$ or $-1$ (depending on the spin of the particle). We consider two cases:
If $T^2 |n\rangle = |n\rangle$, we define two states $|n\rangle_\pm = |n\rangle \pm T |n\rangle$. These states satisfy $T |n\rangle_\pm = \pm |n\rangle_\pm$ which is precisely equation (1) with phase $1$ or $-1$. The rest of the proof follows exactly showing that the wave-function $\psi_{\pm,n}(x)$ is a phase times a real function.
If $T^2 |n\rangle = - |n\rangle$, it is not possible to define any linear combination such that (1) is true. The theorem completely fails in this case.
Let me make a comment about the phase $\theta_n$. This phase freedom is actually not physical. This is because, we can rescale the state by an unphysical phase $|n\rangle \to e^{i\alpha_n} |n\rangle$. In this case, (1) is modified to
$$
T ( e^{i\alpha_n} |n\rangle ) = e^{i\theta_n} ( e^{i\alpha_n} |n\rangle ) \implies e^{-i\alpha_n} T |n\rangle = e^{i( \theta_n + \alpha_n)} |n\rangle \implies T |n\rangle = e^{i( \theta_n + 2\alpha_n)} |n\rangle
$$
We used the fact that $T$ is an anti-linear operator so $T(c |n\rangle) = c^* T |n\rangle$. We can now see that if we choose $\alpha_n = - \frac{1}{2} \theta_n$, we can remove the phase in equation (1). In terms of the new phase rescaled state, the equation is simply $T|n\rangle = |n\rangle$ which then implies that $\psi_n(x)$ is real.
