According to Wikipedia, if a system has $50\%$ chance to be in state $\left|\psi_1\right>$ and $50\%$ to be in state $\left|\psi_2\right>$, then this is a mixed state.

Now, consider the state 
$$\left|\Psi\right>=\frac{\left|\psi_1\right>+\left|\psi_2\right>}{\sqrt{2}},$$ which is a superposition of the states $\left|\psi_1\right>$ and $\left|\psi_2\right>$. Let $\left|\psi_i\right>$ be eigenstates of the Hamiltonian operator. Then measurements of energy will give $50\%$ chance of it being $E_1$ and $50\%$ of being $E_2$. But this then corresponds to the definition above of mixed state! However, superposition is defined to be a pure state.

So, what is the mistake here? What is the real difference between mixed state and superposition of pure states?

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.
