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.
I'm answering only part of your question, please bear in mind that some of the topics that you pointed out are doubts I have myself.
Concerning the comment from Feynman, I usually don't like the authority arguments. If your only reason to believe something is that someone said it, it's not a good reason. Just as some examples, Newton believed that light had no wave-like behaviour, which is simply wrong.
About the variational formulation of Schrödinger Equation. If by variational you mean, having a Lagrangian, the following Lagrangian does the job:
$
\mathcal L = \frac{i\hbar}{2}\left(\psi^*\partial_t \psi - \psi \partial_t \psi^*\right) - \frac{\hbar^2}{2m}\left(\nabla \psi\right)\cdot \left(\nabla \psi^*\right) - V\psi^*\psi
$
With the assumption that you have a complex classical scalar field, $\psi$, and that you can calculate the Euler-Lagrange equation separately for $\psi$ and $\psi^*$. As with anything using the extremum action principle, you really have to guess which your Lagrangian is based in Symmetry principles + some guideline for your problem, here is no different.
I don' t know how did Schrödinger made the original derivation, but if you make a slight change of variables in the above lagrangian, you get very close to what you have written above.
The trick is writing $\psi = \sqrt n e^{iS/\hbar}$, where $n$ is the probability density and $S$ is, essentially the phase of the wave-function. If you have trouble with the derivation, I can help you latter.
Again about Feynman's quote. It's possible to arrive at Schrödinger equation without passing through an arbitrary heuristic procedure, that way is developed by Ballentine's book. You still have to postulate where you live, if you consider that arbitrary or not, it is completely up to you.
The other point that it is possible to consistently and rigorously construct a quantum theory from any classical mechanical theory, using quantization by deformation. That's why I think it's a bit false that "It's not possible to derive it from anything you know. It came out of the mind of Schrödinger", as Feynman, and many other great names, said.
Best Answer
None of these things violate the conservation of energy. They merely change energy from one form (e.g. mass) to another form (e.g. electromagnetic). In all these cases the energy is the same before and after the interaction.
The specific form that you posted is the time independent form. That is based on the Hamiltonian, which is more or less as you say, an expression of the total energy. According to Noether’s theorem the conservation of energy is directly tied to time translation invariance. Therefore, it should not be surprising that the time independent form features energy prominently.