Why is the general solution of Schrodinger's equation a linear combination of the eigenfunctions?

You are starting from the incorrect point. The argument follows by linearity of the equation.
Suppose $\Psi_k(x,t)$ is solution of the time dependent Schr$\ddot{\hbox{o}}$dinger equation: $$ i\hbar \frac{\partial }{\partial t}\Psi_k(x,t)=-\frac{\hbar^2}{2m}\frac{\partial^2\Psi_k(x,t)}{\partial x^2}+U(x)\Psi_k(x,t)\, . $$ Then: $$ \Phi(x,t)=a_1\Psi_1(x,t)+a_2\Psi_2(x,t) $$ is also a solution since $$ i\hbar \frac{\partial }{\partial t}\Phi(x,t) =a_1\left(i\hbar \frac{\partial }{\partial t}\Psi_1(x,t)\right)+a_2 \left(i\hbar \frac{\partial }{\partial t}\Psi_2(x,t)\right) $$ and \begin{align} -\frac{\hbar^2}{2m}\frac{\partial^2\Phi(x,t)}{\partial x^2}+U(x)\Phi(x,t) &=a_1\left(-\frac{\hbar^2}{2m}\frac{\partial^2\Psi_1(x,t)}{\partial x^2}+U(x)\Psi_1(x,t)\right)\\ &\quad + a_2\left(-\frac{\hbar^2}{2m}\frac{\partial^2\Psi_2(x,t)}{\partial x^2}+U(x)\Psi_2(x,t)\right)\, . \end{align} These follow simply from the known rule valid for any two differentiable functions $f$ and $g$: $\partial (f+g)/\partial t=\partial f/\partial t+\partial g/\partial t$, and similarly for the partials w/r to $x$. Combining these last two equations you get an identity for any $a_1$ and $a_2$ since each $\Psi_k(x,t)$ is independently a solution. Of course this simply extends to an arbitrary number of terms in the linear combination.

Note the eigenvalue of the time-independent part never enters in this argument. The final step is to observe that separation of variables in the time-dependent equation yields $\Psi_k(x,t)=e^{-iE_k t}\psi_k(x)$ with $\psi_k(x)$ an eigenfunction of the time-independent equation, but again, this does not enter in the argument.


Edit: note this is in contradistinction with the time-independent equation. When $$ -\frac{\hbar^2}{2m}\frac{d^2\psi_k(x)}{dx^2}+U(x)\psi_k(x)=E_k\psi_k(x) $$ the the right hand side is must a multiple of the original function. With this observation, note then that a linear combination $$ \psi(x)=a_1\psi_1(x)+a_2\psi_2(x) $$ will in general NOT be a solution of the time-independent equation because \begin{align} \left(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+U(x)\right)\psi(x) &=a_1\left(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+U(x)\right)\psi_1(x)\\ &\qquad+a_2\left(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+U(x)\right)\psi_2(x) \\ &=a_1E_1\psi_1(x)+a_2E_2\psi_2(x)\\ &=E_1(a_1\psi_1(x)+a_2\psi_2(x))+(E_2-E_1)a_2\psi_2(x)\\ &=E_1\psi(x)+(E_2-E_1)a_2\psi_2(x) \end{align} will NOT be a multiple of $\psi(x)$ unless $E_1=E_2$.


Maybe to answer your question it is useful to start from a slightly different perspective, conceptually. In Quantum Mechanics a system is described by giving a Hilbert space, whose vectors represent states of the system (actually these are only part of the "states" called pure states, and any two vectors that are complex multiples of each other represent the same state. Don't worry about this for the moment) and a prescription of how the state of system evolves with time. In this case this prescription is the Schroedinger equation, which I write as an equation of Hilbert space vectors (in your case, the Hilbert space is the space of square integrable functions over $\mathbb{R}$, called $L^2(\mathbb{R})$) $$i\hbar\partial_t\psi_t=\hat{H}[\psi_t]$$

The operator $\hat{H}$ is self adjoint and mathematics tells us that in this case the unique solution to this equation with given initial condition $\psi_0$ (i.e. the state your system starts in) is $$\psi_t=e^{-\frac{i}{\hbar}\hat{H}t}[\psi_0]$$

where $U_t:=e^{-\frac{i}{\hbar}\hat{H}t}$ is a unitary operator called "time evolution" (compare this to the fact that the matrix exponential of $iA$, where $A$ is a hermitian matrix, is a unitary matrix).

Now the entire remaining work is to calculate $U_t\psi_0$ for given $\psi_0$. This is unfortunately really hard in most cases! That's where the time-independent Schroedinger equation comes in. Suppose you have a bunch of eigenvectors $\phi_1, \phi_2 \dots$ of $\hat{H}$, that is they satisfy $H\phi_i=E_i\phi_i$ (I hope there is no confusion when subscripts designate time and when they designate which $\phi$ I'm talking about). Again, mathematics tells us (and it's easy to verify this for matrices, that is for finite dimensional Hilbert spaces) that $$U_t[\phi_i]=e^{-\frac{i}{\hbar}\hat{H}t}[\phi_i]=e^{-\frac{i}{\hbar}E_it}\phi_i$$

This is really nice, because now we don't have to calculate the exponential of an operator. Also all these operators are linear, so if we manage to write our initial state $\psi_0$ as some linear combination of eigenvectors $$\psi_0=\sum_{i=1}^{?} c_i\phi_i$$ for some constants $c_i$ the desired solution takes the form

$$\psi_t=\sum_{i=1}^{?} c_ie^{-\frac{i}{\hbar}E_it}\phi_i$$ Thus we have solved the problem if we can find a way to express any initial condition we might want as a linear combination of eigenvectors of $\hat{H}$. We want to find an orthonormal basis of the Hilbert space consisting of such eigenvectors, then we can express ANY vector as an infinite linear combination.