The state
\begin{equation}
|\Psi \rangle = \frac{1}{\sqrt{2}}\left(|\psi_1\rangle +|\psi_2\rangle \right)
\end{equation}
is a pure state. Meaning, there's not a 50% chance the system is in the state $|\psi_1\rangle$ and a 50% it is in the state $|\psi_2\rangle$. There is a 0% chance that the system is in either of those states, and a 100% chance the system is in the state $|\Psi\rangle$.
The point is that these statements are all made before I make any measurements.
It is true that if I measure the observable corresponding to $\psi$ ($\psi$-gular momentum :)), then there is a 50% chance after collapse the system will end up in the state $|\psi_1\rangle$.
However, let's say I choose to measure a different observable. Let's say the observable is called $\phi$, and let's say that $\phi$ and $\psi$ are incompatible observables in the sense that as operators $[\hat{\psi},\hat{\phi}]\neq0$. (I realize I'm using $\psi$ in a sense you didn't originally intend but hopefully you know what I mean). The incompatibliity means that $|\psi_1 \rangle$ is not just proportional to $|\phi_1\rangle$, it is a superposition of $|\phi_1\rangle$ and $|\phi_2\rangle$ (the two operators are not simulatenously diagonalized).
Then we want to re-express $|\Psi\rangle$ in the $\phi$ basis. Let's say that we find
\begin{equation}
|\Psi\rangle = |\phi_1\rangle
\end{equation}
For example, this would happen if
\begin{equation}
|\psi_1\rangle = \frac{1}{\sqrt{2}}(|\phi_1\rangle+|\phi_2\rangle)
\end{equation}
\begin{equation}
|\psi_2\rangle = \frac{1}{\sqrt{2}}(|\phi_1\rangle-|\phi_2\rangle)
\end{equation}
Then I can ask for the probability of measuring $\phi$ and having the system collapse to the state $|\phi_1\rangle$, given that the state is $|\Psi\rangle$, it's 100%. So I have predictions for the two experiments, one measuring $\psi$ and the other $\phi$, given knowledge that the state is $\Psi$.
But now let's say that there's a 50% chance that the system is in the pure state $|\psi_1\rangle$, and a 50% chance the system is in the pure state $|\psi_2\rangle$. Not a superposition, a genuine uncertainty as to what the state of the system is. If the state is $|\psi_1 \rangle$, then there is a 50% chance that measuring $\phi$ will collapse the system into the state $|\phi_1\rangle$. Meanwhile, if the state is $|\psi_2\rangle$, I get a 50% chance of finding the system in $|\phi_1\rangle$ after measuring. So the probability of measuring the system in the state $|\phi_1\rangle$ after measuring $\phi$, is (50% being in $\psi_1$)(50% measuring $\phi_1$) + (50% being in $\psi_2$)(50% measuring $\phi_1$)=50%. This is different than the pure state case.
So the difference between a 'density matrix' type uncertainty and a 'quantum superposition' of a pure state lies in the ability of quantum amplitudes to interfere, which you can measure by preparing many copies of the same state and then measuring incompatible observables.
Almost. You do know its state. It's $(c_1\Psi_1+c_2\Psi_2)$ (apart from a normalization constant). Remember that your choice of basis vectors to represent the degenerate subspace is arbitrary. There's nothing in the physics distinguishing the two you happened to choose. That superposition state is just as good as any of the other states in that degenerate subspace. This is what it means to project onto the degenerate subspace and renormalize.
Now here's a messier variant on the problem: Suppose the spectrum of the Hamiltonian is continuous, and the energy measurement has a finite precision (with a known error distribution). Now what's the state after the measurement?
There actually is a good answer to this, but it's unlikely to be taught in most courses. You basically end up with a quantum mechanical version of Bayes' theorem. The basic version you're taught in a course is just the limit of an ideal measurement, such that the state of the measurement device is guaranteed to 100% correlate with the thing being measured.
Best Answer
A slight expansion on danimal's comment: you can generally get the state $\psi(x,t)$ from the $\psi(x,0)$ you provided by operating on it with the unitary time evolution operator $\exp(-i \hat{H} t/\hbar)$. Since you know the eigenstates, you can write the Hamiltonian in a diagonal basis and this operator will appear to multiply $\psi_n$ simply by $e^{-i E_n t/\hbar}$.
To find the expectation value of any Hermitian observable $\hat{A}$ corresponding to your state at time $t$, you can simply compute $A(x,t) = \langle \psi(x,t) \vert \hat{A} \vert \psi(x,t)\rangle$. To find the time average, you could simply treat this expectation value as a classical function and average it over an interval, say by integration.