For Hamiltonian operators of the form
$$H = -\frac{\hbar^2}{2m} \Delta + V(x)$$
$V$ must be real to assure that $H$ is at least symmetric and eigenvalues $E$ are real. Therefore, if $0\neq \psi= \psi(x)$ is an eigenfunction, so that:
$$-\frac{\hbar^2}{2m} \Delta \psi + V(x) \psi = E\psi \tag{1}$$
taking the complex conjugate of both sides you have:
$$-\frac{\hbar^2}{2m} \Delta \overline{\psi} + V(x)\overline{\psi} = E\overline{\psi}\:.\tag{2}$$
(1) and (2) together imply that the real valued functions $\Psi= \psi + \overline{\psi}$ and $\Phi= i(\psi - \overline{\psi})$ are eigenfunctions with the same eigenvalue $E$. Notice that the pair $\Psi$, $\Phi$ encompass the same information as that of the pair $\psi$ and $\overline{\psi}$: They are linearly independent if and only if $\psi$ and $\overline{\psi}$ are and generate the same vector subspace (of the eigenspace associated to $E$).
If an eigenspace with energy $E$ has dimension $\geq 2$, you may have a
nonvanishing flux of a generic eigenfunction with that energy.