There is no need for the solution $\psi(x)$ to be real. What must be real is the probability density that is "carried" by $\psi(x)$. In some loose and imprecise intuitive way, you may think about a TV image carried by electromagnetic waves. The signal that travels is not itself the image, but it carries it, and you can recover the image by decoding the signal properly.
Somewhat similarly, the complex wave function that is found by solving Schrödinger equation carries the information of "where the particle is likely to be", but in an indirect manner. The information on the probability density $P(x)$ of finding the particle is recovered from $\psi(x)$ simply by multiplying it times its complex conjugate:
$$\psi(x)^*\psi(x) = P(x)$$
that gives a real function as a result. Note that it is a density: what you compute eventually is the probability of finding the particle between $x=a$ and $x=b$ as $\int_{a}^{b} P(x) dx$
As you know, when you multiply a complex number(/function) times its complex conjugate, the information on the phase is lost:
$$\rho e^{i \theta}\rho e^{-i \theta}=\rho^{2}$$
For that reason, in some places one can (not quite correctly) read that the phase has no physical meaning (see footnote), and then you may wonder "if I eventually get real numbers, why did not they invent a theory that directly handles real functions?".
The answer is that, among other reasons, complex wave functions make life interesting because, since the Schrödinger equation is linear, the superposition principle holds for its solutions. Wave functions add, and it is in that addition where the relative phases play the most important role.
The archetypical case happens in the double slit experiment. If $\psi_{1}$ and $\psi_{2}$ are the wave functions that represent the particle coming from the hole number $1$ and $2$ respectively, the final wave function is
$$\psi_{1}+\psi_{2}$$
and thus the probability density of finding the particle after it has crossed the screen with two holes is found from
$$P_{1+2}= (\psi_{1}+\psi_{2})^{*}(\psi_{1}+\psi_{2}) $$
That is, you shall first add the wave functions representing the individual holes to have the combined complex wave function, and then compute the probability density. In that addition, the phase informations carried by $\psi_{1}$ and $\psi_{2}$ play the most important role, since they give rise to interference patterns.
Comment: Feynman is quoted to have said "One of the miseries of life is that everybody names things a little bit wrong, and so it makes everything a little harder to understand in the world than it would be if it were named differently." It is quite similar here. Every book says that the phase of the wave function has no physical meaning. That is not 100% correct, as you see.
First, since the potential is unbounded from below at the left side, the particle has continuous spectrum. This means that what you just have to do is compute the coefficients $c_i$, taking $E$ (real, never complex!) as input parameter.
One of these coefficients, $c_3$, is in fact arbitrary, because it only influences normalization, not smoothness of wavefunction. Thus, to compute the connection coefficients you have to find $c_1,c_2,c_4,c_5$.
And finally, your equation is a simple system of $4$ algebraic equations with $4$ unknowns. If I take $c_3=1$ and use arbitrary letters to denote all those Airy functions and their derivatives, I get
$$\left\{\begin{align}
c_4A+c_5B&=c_1C+c_2D,\\
c_4G+c_5H&=c_1I+c_2J,\\
c_1K+c_2L&=M,\\
c_1N+c_2O&=P.
\end{align}\right.$$
Solve it for $c_i$ and you've solved almost the whole of your problem. Now if you need normalized wavefunctions, you have to use some scheme for normalization of unbounded states, e.g. so called "normalization by Dirac delta", discussed e.g. in [1].
References:
- L.D. Landau & E.M. Lifshitz, Quantum Mechanics. Non-Relativistic Theory, $\S5$.
Best Answer
In ordinary quantum mechanics, two wavefunctions represent the same physical state if and only if they are multiples of each other, that is $\psi$ and $c\psi$ represent the same state for any $c\in\mathbb{C}$. If you insist on wavefunctions being normalized, then $c$ is restricted to complex numbers of absolute value 1, i.e. number of the form $\mathrm{e}^{\mathrm{i}k}$ for some $k\in[0,2\pi)$.
If $\psi(x) = a(x) + \mathrm{i}b(x)$ is a solution of the Schrödinger equation, then it is not automatically true that $a(x)$ and $b(x)$ are also solutions. It is "accidentally" true for the time-independent Schrödinger equation because applying complex conjugation shows us directly that $\psi^\ast(x)$ is a solution if $\psi(x)$ is, and $a(x)$ and $b(x)$ can be obtained by linear combinations of $\psi(x)$ and $\psi^\ast(x)$.
There are now two cases: If $\psi$ and $\psi^\ast$ are not linearly independent - i.e. one can be obtained from the other by multiplication with a complex constant - then the space of solutions for this energy is still one-dimensional, and there's only a single physical state. If they are linearly independent, then there are at least two distinct physical states with this energy.
Note that already the free particle with $V=0$ gives a counter-example to the claim that all solutions for the same energy have the same values for all expectation values. There we have plane wave solutions $\psi(x) = \mathrm{e}^{\mathrm{i}px}$ and $\psi^\ast(x) = \mathrm{e}^{-\mathrm{i}px}$ that are linearly-independent complex conjugates with the same energy that differ in the sign of their expectation value for the momentum operator $p$.