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.
Start with the time-dependent Schrödinger equation
$$ \hat H \Psi(r, t) = i \hbar \partial_t \Psi (r, t). $$
Using the ansatz $\Psi(r,t) = \psi(r) f(t)$ yields
$$ f(t) \hat H \psi(r) = i \hbar \psi(r) \partial_t f(t)$$
and, via the standard separation of variables trick,
$$ i\hbar\frac{\dot f(t)}{f(t)} = \text{const} = \frac{\hat H\psi(r)}{\psi(r)}; $$
the two sides are equal, but the LHS depends only on $t$ while the RHS depends only on $r$, so they each have to be constant. Let us call this constant $E$. Then, for the time-dependent part, we get
$$ \dot f(t) = -i \frac{E}{\hbar} f(t) $$
Which is manifestly solved by
$$ f(t) = \exp\left\{-i \frac{E}{\hbar} t\right\}. $$
For the spatial part, we find the time-independent Schrödinger equation
$$ \hat H \psi(r) = E \psi(r), $$
which as you observe can be viewed as an eigenvalue equation for $\hat H$. This motivates the choice of name for $E$: Physically, the eigenvalue of the Hamiltonian is the energy.
Best Answer
No.
If we view the Schrödinger equation simply as a second order differential equation, then no, it doesn't have a general solution that we can find. Let me add that, from the mathematical standpoint, the equation alone is not even sufficient to have a well posed problem; it needs to be supplemented by the boundary conditions. The solutions for the same type of potential may be easy to get for one set of boundary conditions and hard or even impossible for the other. For example, the harmonic potential case, $V(x) = \frac{m\omega^2x^2}{2}$ is exactly solvable with open boundary conditions, but solving it with hard wall conditions, e.g., at $x= \pm a$, already requires solving transcendental equations; i.e. the solution is no longer exact. (However, some particular cases, such as the case of one hard wall at $x=0$ are still doable.)
But yes.
One may question the very notion of an exact solution: usually it means a solution in terms of simple functions that can be handled with pan and paper. Some people would generalize this to using special functions — Bessel functions, hypergeometric functions, etc. However, this notion of "exact" ultimately traces itself to our ability to calculate the numbers. When we need to calculate a sine, an exponent, or a Bessel function, most of us turn to a computer... but once we accept that computers can be used, almost any one particle Schrödinger equation is solvable. Note that this is not true for many-particle problems, where the problems easily get NP-hard, i.e. impossible using modern computational power (but possibly solvable with a quantum computer).
And formally, yes.
Using the evolution operator one can construct a formally exact solution, evolving from an exactly solvable case, which is quite a common procedure when developing all kinds of perturbation expansions.