You should be careful not to confuse probabilities and probability densities. Let's start with the wavefunction corresponding to your example:
$$\psi(x)=\sqrt{\tfrac 1 3}\psi_0(x)+\sqrt{\tfrac 2 3}\psi_1(x)$$
We can now ask two kinds of questions (we can ask many more questions, but let's focus on these two):
- When I measure the energy of this particle, what is the probability of getting $E_0$ or $E_1$?
- What is the probability (density) of finding the particle at a certain location?
Probability of finding $E_0$ or $E_1$
To answers question 1 we can make use of the orthogonality of the wavefunctions. The inner product between two wavefunctions is defined as
$$\langle \psi|\phi\rangle=\int\mathrm d x\, \psi^*(x)\phi(x).$$
The energy eigenstates form an orthonormal basis$^\dagger$, so
\begin{align}
\langle \psi_0|\psi_0\rangle=1\\
\langle \psi_1|\psi_1\rangle=1\\
\langle \psi_0|\psi_1\rangle=0
\end{align}
The probability of measuring energy $E_0$ is now given by
\begin{align}
p(E_0) &= \left|\langle \psi_0|\psi\rangle\right|^2\\
&=\left|\sqrt{\tfrac 1 3}\underbrace{\langle \psi_0|\psi_0\rangle}_{1} + \sqrt{\tfrac 2 3}\underbrace{ \langle \psi_0|\psi_1\rangle }_{0}\right|^2\\
&=\frac 1 3
\end{align}
Similarly, $p(E_1)=\frac 2 3$.
Probability density of $\psi_0 + \psi_1$
Now what does the expression $|\psi(x)|^2$ mean before we integrate? The interpretation is that $|\psi(x)|^2$ tells you how likely it is to find a particle at $x$. Well, the probability of finding it at $x$ is always zero, but if you integrate over a small region you get
$$p(a\leq x\leq b)=\int\mathrm d x\,|\psi(x)|^2.$$
If you are confused as to why the probability at location $x$ is zero, you are not alone. Search on google for 'probability density' to find nice explanations.
When you add two eigenstates together, the probability density functions are modified. This is called interference. The probability density of the superposition of two functions is not simply the sum of the probability densities. As you showed, $|\psi|^2$ is given by
$$|\psi|^2=\tfrac 1 3|\psi_0|^2+\tfrac 2 3|\psi_1|^2+\tfrac {2\sqrt 2} 3\text{Re}\left[\psi_0^*(x)\psi_1(x)\right]$$
The terms $|\psi_0|^2$ and $|\psi_1|^2$ are just the probabilities densities added together, as you would expect naively. The term involving $\psi_0\psi_1$ is an interference term. It can be negative and it is the root cause of interference patterns like in the double slit experiment.
Here is an example of how that would look like.
$\dagger$ The reason that these states are orthonormal is because they are the eigenstates of a Hermitian operator. To see this for yourself, you can create a random matrix $M$ in Mathematica/Matlab/Python, then create a Hermitian matrix using $H=M+M^\dagger$. If you calculate the eigenvectors, they will all be orthogonal with respect to each other. They might even be normalized if your solver is nice.
