I understand that the Schrodinger equation is actually a principle that cannot be proven. But can someone give a plausible foundation for it and give it some physical meaning/interpretation. I guess I'm searching for some intuitive solace here.
[Physics] What does the Schrodinger Equation really mean
complex numbersquantum mechanicsschroedinger equationwave-particle-dualitywavefunction
Related Solutions
So by "general solution" you mean something like $$ \psi_{general} = \sum_{n = 1}^\infty c_n \left[\sqrt{\frac{2}{a}}\sin\left(\frac{n\pi x}{a}\right)\right],\ \ c_n \in \mathbb{C}, $$ right? (I'm supposing that the infinite square well's walls are at $0$ and $a$, and I'm writing $c_j$ for arbitrary complex numbers. Also, note that I've pulled a constant factor $\sqrt{2/a}$ out of the constants and put it into brackets---you'll see the reason for that in a bit.)
The thing is, this expression $\psi_{general}$ doesn't have any physical meaning as is. You also have to pick values for the $c_n$ before you can think of it as a particular state. For example you could write $$ \psi = \sum_{n = 1}^\infty c_n\left[\sqrt{\frac{2}{a}}\sin\left(\frac{n\pi x}{a}\right)\right],\\ c_n = 1/\sqrt{100},\ \ n \le 100, \text{$n$ odd}\\ c_n = 0, \ \ n \le 100, \text{$n$ even}\\ c_n = 0, \ \ n \ge 100 $$ or just $$ \psi = \frac{1}{\sqrt{100}}\sum_{n = 1}^{50} \left[\sqrt{\frac{2}{a}}\sin\left(\frac{(2n - 1 )\pi x}{a}\right)\right]. $$
Vladimir Kalitviansky is exactly right that this sort of state is a "superposition" of a whole bunch of other states. What that means, though, is pretty deep. The idea is that when you're in a superposition of many states, you have a certain probability to end up in any of them. If I measure the energy of a particle in state $\psi$, the wavefunction will collapse into a randomly selected one of the states $$ \psi_n = \sqrt{\frac{2}{a}}\sin\left(\frac{n\pi x}{a}\right) $$ (I think these are the separable solutions to which you refer) and you'll measure an energy $E_n = \frac{n^2\pi^2\hbar^2}{2ma^2}$; the probability that the state will end up in any particular state $n$ is $|c_n|^2$, or in this case $1/100$. So: you start out in $\psi$, I measure your energy to be $E_n$, and you end up in a state $\psi_n$.
It's important to understand that which state you end up in is totally random. There is absolutely no way to tell beforehand what state it'll be, and therefore what energy will be measured---no hidden variables that determine what the state will turn out, and no correlations with other variables. Sure, you know that some states (those with the largest $c_n$ associated) will be more likely than others, but that's it. This is the principle behind many hardware random number generators.
You're presumably wondering about those functions $\psi_n$. I haven't told you why the state will collapse into those particular functions. What's so special about the $\psi_n$ as opposed to some other function---say, $\arccos(\sinh(\exp(x)))$? Also, I just pulled the numbers $E_n$ out of nowhere, which is distinctly unsporting.
What's going on is that these functions are eigenfunctions of the Hamiltonian with eigenvalue $\frac{n^2\pi^2\hbar^2}{2m}$. The vocabulary about eigen-stuff is just putting names to some facts I you've already realized: namely, that when I apply the Hamiltonian $$ H = -\frac{\hbar^2}{2m} \frac{d^2}{dx} $$ to one of the states $\psi_n$, I get back $\psi_n$ (the eigenfunction) multiplied by a constant $E_n$ (the eigenvalue)---that is to say, $$ -\frac{\hbar^2}{2m} \frac{d^2}{dx} \psi_n = \frac{n^2\pi^2\hbar^2}{2ma^2}\sqrt{\frac{2}{a}} \sin\left(\frac{n\pi x}{a}\right) = E_n \psi_n, $$ which is just the (time-independent) Schrodinger equation. So, t'm in ao rephrase my summary: you start out in a superposition of states, I measure your energy, and you end up in some eigenfunction of the Hamiltonian.
This is one (well, two) of the general axioms of quantum mechanics: if the eigenfunctions of the Hamiltonian are $\psi_n$ (these don't have to be the $\psi_n$ I wrote down above, if the Hamiltonian is not the square-well Hamiltonian) and I measure the energy of a particle in a state $\psi$, then the particle will end up in a state $\psi_n$ with probability $$ P(n) = \left|\int_{-\infty}^\infty \psi_n^* \psi\ dx\right|^2 = |c_n|^2 $$ (where $^*$ denotes complex conjugation---that always causes me some confusion, since mathematicians use the symbol differently). That last step (where I bring in the $c_n$s) is not totally obvious. It follows from the fact that the $\psi_n$ are orthonormal---that is, $\int_{-\infty}^\infty \psi_j^*\psi_k\ dx = 1$ if $j = k$, and $0$ if $j \ne k$. You'll want to work out for yourself exactly how that works. First, notice that the $\psi_{general}$ can be written nicely in terms of the $\psi_n$; then, plug that $\psi_{general}$ into the above expression for $P(n)$; finally, prove that the particular $\psi_n$ I wrote for the square-well case are orthonormal using a trig identity (that's why I put the constant out front of them).
This imposes the very important constraint that $$ \sum |c_n|^2 = \sum P(n)= 1. $$ The particle has to end up in at least one of the states, so the sum of the probabilities can't be less than $1$, and it would be total nonsense for the total probability to be greater than $1$, so $1$ it is. This is called "normalization", as in "the state is normalized".
As another fun exercise, explore the state I gave you as a specific example. Plot it ($x$ in units of $a$), change that 100 to a bunch of different numbers $N$ (you'll have to change the value of the $c_j$ so that the state is normalized), and plot again. You should see that the state is more or less localized around $x = a/2$. What is the width of the state, and how does that width (call it $\Delta x$) change with N? Given that the momentum of an eigenstate $\psi_n$ is $p_n = \sqrt{2mE} = n\pi\hbar /a$, how does this relate to Heisenberg's uncertainty principle?
I should note that I've skipped a whole bunch of details, and not just in that exercise at the end. If any of them are bothering you, ask! Hope this helps.
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
This is a fairly basic approach suitable for students who have finished at least one semester of introductory Newtonian mechanics, are familiar with waves (including the complex exponential representation) and have heard of the Hamiltonian at a level where $H = T + V$. As far as I understand it has no relationship to Schrödinger's historical approach.
Let's take up Debye's challenge to find the wave equation that goes with de Broglie waves (restricting ourselves to one dimension merely for clarity).
Because we're looking for a wave equation we will suppose that the solutions have the form $$ \Psi(x,t) = e^{i(kx - \omega t)} \;, \tag{1}$$ and because this is suppose to be for de Broglie waves we shall require that \begin{align} E &= hf = \hbar \omega \tag{2}\\ p &= h\lambda = \hbar k \;. \tag{3} \end{align}
Now it is a interesting observation that we can get the angular frequency $\omega$ from (1) with a time derivative and likewise wave number $k$ with a spacial derivative. If we simply define the operators1 \begin{align} \hat{E} = -\frac{\hbar}{i} \frac{\partial}{\partial t} \tag{4}\\ \hat{p} = \frac{\hbar}{i} \frac{\partial}{\partial x} \; \tag{5}\\ \end{align} so that $\hat{E} \Psi = E \Psi$ and $\hat{p} \Psi = p \Psi$.
Now, the Hamiltonian for a particle of mass $m$ moving in a fixed potential field $V(x)$ is $H = \frac{p^2}{2m} + V(x)$, and because this situation has no explicit dependence on time we can identify the Hamiltonian with the total energy of the system $H = E$. Expanding that identity in terms of the operators above (and applying it to the wave function, because operators have to act on something) we get \begin{align} \hat{H} \Psi(x,t) &= \hat{E} \Psi(x,t) \\ \left[ \frac{\hat{p}^2}{2m} + V(x) \right] \Psi(x,t) &= \hat{E} \Psi(x,t) \\ \left[ \frac{1}{2m} \left( \frac{\hbar}{i} \frac{\partial}{\partial x}\right)^2+ V(x) \right] \Psi(x,t) &= -\frac{\hbar}{i} \frac{\partial}{\partial t} \Psi(x,t) \\ \left[ -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2}+ V(x) \right] \Psi(x,t) &= i\hbar \frac{\partial}{\partial t} \Psi(x,t) \;. \tag{6}\\ \end{align} You will recognize (6) as the time-dependent Schrödinger equation in one dimension.
So the motivation here is
but this is not anything like a proof because the pass from variable to operators is pulled out of a hat.
As an added bonus if you use the square of the relativistic Hamiltonian for a free particle $(pc)^2 - (mc^2)^2 = E^2$ this method leads naturally to the Klein-Gordon equation as well.
1 In very rough language an operator is a function-like mathematical object that takes a function as an argument and returns another function. Partial derivatives obviously qualify on this front, but so do simple multiplicative factors: because multiplying a function by some factor returns another function.
We follow a common notational convention in denoting objects that need to be understood as operators with a hat, but leaving the hat off of explcit forms.