Solving the 1d Schroedinger Equation numerically. How to get complex phases? There is a nice method that lets you solve the Schroedinger Equation in 1d numerically, by discretizing space. Tranforming the Eigenfunction problem into an eigenvector problem with a finite-dimensional matrix that can be diagonalized. After finding the eigenvalues (Eigenenergies), the eigenvectors can be computed directly. The i-th entry of one of these eigenvectors is the amplitude of the corresponding wavefunction at position $x_i$.
The algorithm is explained here by a stackexchange user (he wrote this in response to a question regarding a similar topic):

I'm glad I can finally answer this.
Numerov's method as described on Wikipedia is not how you want to
proceed here. To give you an idea of how to proceed, let's start with
a simplified version of the method. What I'm going to do is to just
naively discretize the differential operator, like so:
$$ \frac{d^2}{dx^2}\psi \approx \frac{\psi(x-d)-2\psi(x)
  +\psi(x+d)}{d^2} \equiv \frac{\psi_{n-1} - 2 \psi_n + \psi_{n+1}}{d^2}$$
The last equation is just a definition -- I'm treating space as if
it's a lattice with points $d$ apart and I'm calling $\psi_n$ the
value of the wavefunction on the $n$th point.. Now the Schrödinger
equation reads in this notation:
$$-\frac{\hbar}{2m}\frac{\psi_{n-1} - 2 \psi_n + \psi_{n+1}}{d^2} +
  V_n \psi_n = E \psi_n$$
But this is a matrix equation! Let me be more explicit:
$$\begin{align}
  -\frac{\hbar}{2md^2}&\left(\begin{array}{cccccccc}
  -2 &  1   &      &      &    \\ 1  & -2   & 1    &      &    \\    &\ddots&\ddots&\ddots&    \\    &      & 1    & -2   & 1  \\    &     
  &      &  1   & -2 \\ \end{array}\right) \left(\begin{array}{c} \psi_1
  \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) +\\
  &\qquad\qquad\quad\left(\begin{array}{cccccccc} V_1 &     &         & 
  &    \\
      & V_2 &         &         &    \\
      &     &  \ddots &         &    \\
      &     &         & V_{N-1} &    \\
      &     &         &         & V_N\\ \end{array}\right)\left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\
  \psi_{N-1} \\ \psi_N \end{array}\right) = E \left(\begin{array}{c}
  \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right)
  \end{align}$$
Of course, I had to pick some integer $N$ that just corresponds to
placing the system in a box that's "big enough" in order to have a
finite system.
It is clear now that what you have is a matrix eigenvalue problem of
the form
$$A\psi = E \psi$$
and you may proceed to diagonalize it in whatever way you choose. Note
that we call $A$ a tridiagonal matrix, for obvious reasons. You might
want to take care of the boundary conditions first -- you do that by
setting $\psi_1 = \psi_n = 0$ before you construct the matrix, which
corresponds to setting the first and last columns to zero. This will
give you some spurious eigenfunctions with zero eigenvalue that you
can just discard. If you have different boundary conditions, you're
out of luck -- I don't know of a way that makes it work in general.
Now you just have to redo the above with the full Numerov's method,
which will be a bit more complicated, and you're all set.
This seems to be the way you're looking for but of course it's not the
only way to do this. Griffiths describes one called "wag the dog" that
consists of the following: you guess an initial value for an energy
and compute the wavefunction however you want (RK, for instance).
Chances are it won't be an eigenvalue so the wavefunction will blow up
at infinity. It might go to $+\infty$ for large $x$ or it might go to
$-\infty$. You now slowly vary the energy until the behavior at
infinity "flips", that is, until the tail "wags". That will allow you
to constrain the value of a single energy eigenvalue and give you the
form of the wavefunction out to some maximum size that increases as
the chosen $E$ approaches the exact eigenvalue.

For a real physical potential, the matrix is symmetric and real. As far as I'm aware this means that the eigenvectors have all real entries (or all complex, but then one can extract a global phase). Analytical solutions to certain potential however do have complex wavefunctions.
The procedure only works with the following inital condition: The eigenfunctions need to be zero at the boundaries (e.g. in the setting of an "infinite well")
Is this the reason why the wavefunctions are all real?
When looking at practical problems this should not be much of an issue as one could just define the relevant potential in the middle far away from the boundaries.
 A: For Hamiltonians like these with just a position dependent potential,
taking the complex conjugate of the Schrodinger differential equation
shows that the complex conjugate is also an eigenstate with the same
energy. Therefore the real part and the imaginary parts are also
eigenstates. Therefore, for these problems, you can always choose the
energy eigenstates to be real.
Many physical systems correspond to  particles in a wave packet with
some momentum or angular momentum. For example for a free particle you
likely want to expand in plane waves $e^{ikx}$ and $e^{-ikx}$
rather than $\cos(kx)$ and $\sin(kx)$. So you choose these complex solutions.
Since momentum is not conserved with your hard wall boundary conditions
to get these solutions, you would need to make your system periodic.
That is the first and last rows need to be replaced with
$-2\ 1 \ ...0\ 1$ and $1\ 0 ... 1 \ -2$. You would then get degenerate
solutions. Your computer program likely will return real eigenvectors,
but you could take two linear combinations of the doubly degenertate
ones to give two independent complex eigenvectors.
