When Eigenfunctions/Wavefunctions are real?
All bound states can typically be chosen to have real-valued wavefunctions. The reason for this is that their wavefunction obeys a real differential equation, $$ -\frac{\hbar^2}{2m}\nabla^2\psi(\mathbf r)+V(\mathbf r)\psi(\mathbf r)=E\psi(\mathbf r)$$ and therefore for any solution you can construct a second solution by taking the complex conjugate $\psi(\mathbf r)^\ast$. This second solution will either be
- linearly dependent on $\psi$, in which case $\psi$ differs from a real-valued function by a phase, or
- linearly independent, in which case you can "rotate" this basis into the two independent real-valued solutions $\operatorname{Re}(\psi)$ and $\operatorname{Im}(\psi)$.
For continuum states this also applies, but things are not quite as clear as the boundary conditions are not invariant under conjugation: incoming scattering waves with asymptotic momentum $\mathbf p$, for example, behave asymptotically as $e^{i\mathbf p\cdot \mathbf r/\hbar}$, and this changes into outgoing waves upon conjugation. Thus, while you can still form two real-valued solutions, they will be standing waves and the physics will be quite different.
In the second case, when you have a degeneracy, the physical characteristics of the real-valued functions are in general different to those of the complex-valued ones. For example, in molecular physics, $\Pi$ states typically have such a degeneracy: you can choose
- $\Pi_x$ and $\Pi_y$ states, which are real-valued, have a node on the $x$ and $y$ plane, resp., have a corresponding factor of $x$ and $y$ in the wavefunction, and have zero expected angular momentum component along the $z$ axis, or
- $\Pi_\pm$ states, which have a complex factor of $x\pm i y$ and no node, and have definite angular momentum of $\pm \hbar$ about the $z$ axis.
Thus: you can always choose a real-valued eigenstate, but it may not always be the one you want.
In addition to Emilio's great answer, and in answer to your second question: Specifically in 1D potential problems (i.e. $\hat H = \frac{1}{2m} \hat p^2 + V(\hat x)$), all the bound states can simultaneously be made real. This is because of the theorem that bound states in 1D are non-degenerate; then, $\psi$ and $\psi^*$, which are both solutions in any dimensionality, must be linearly dependent.
The situation is different if you have a magnetic field; then the prescription ("minimal coupling") is to replace $\hat p \gets (\hat p - e A(\hat x))$, which results in a complex Hamiltonian ($A$ is the vector potential).
Also, the argument above holds for a "well-behaved" potential; see http://arxiv.org/abs/0706.1135 for a modern take on that.
Finally, regarding your comment about a band-structure calculation: I am not sure this is what you mean, but, at least in the context of certain codes, you will often hear it said that a calculation for a centrosymmetric structure without spin-orbit coupling is "real", but this does not mean that the true eigenfunctions $\psi_{nk}$ are real-valued (remember, they contain a plane-wave factor $e^{ikr}$). Rather, the coefficients in the basis used for the calculation can be chosen to be real.