Let's suppose you have two point charges $q_1$ and $q_2$. Each charge has its own electric field $\vec{E}_1$ and $\vec{E}_2$, and the total electric field is $\vec{E} = \vec{E}_1 + \vec{E}_2$. If we look at the total potential energy "stored in the fields" in this configuration, we can split it up into three pieces:
\begin{align*}
U = \iiint u \, dV &= \frac{\epsilon_0}{2} \iiint \left( \vec{E}_1 + \vec{E}_2 \right)^2 \, dV \\
&= \frac{\epsilon_0}{2} \iiint \left( \vec{E}^2_1 + \vec{E}^2_2 + 2 \vec{E}_1 \cdot \vec{E}_2 \right) \, dV \\
&= \underbrace{\frac{\epsilon_0}{2} \iiint \vec{E}^2_1 \, dV}_{U_1} + \underbrace{\frac{\epsilon_0}{2} \iiint \vec{E}^2_2 \, dV}_{U_2} + \underbrace{ \epsilon_0 \iiint \vec{E}_1 \cdot \vec{E}_2 \, dV}_{U_\text{int}}
\end{align*}
Now, what do each of these three pieces mean? $U_\text{int}$ turns out to be the easiest to interpret: if you calculate this integral over all over space assuming that the charges are a distance $r$ apart, you get
$$
U_\text{int} = \frac{q_1 q_2}{4 \pi\epsilon_0 r};
$$
in other words, this is the potential energy that we know & love. (Doing this exercise is a fun way to test your calculational mettle.) This term can therefore be thought of as the energy due to the interaction between the two charges $q_1$ and $q_2$. But what about $U_1$ and $U_2$? It's not too hard to see that these two quantities are, in fact, infinite; writing out the integral in spherical coordinates, we get
$$
U_1 = \frac{\epsilon_0}{2} \iiint_\text{all space} \left( \frac{q_1}{4 \pi \epsilon_0 r^2} \right)^2 r^2 \sin \theta \, dr \, d\theta \, d \phi = \frac{q_1^2}{8 \pi \epsilon_0} \int_0^\infty \frac{1}{r^2} dr = \frac{q_1^2}{8 \pi \epsilon_0} \left[ \frac{1}{r} \right]_0^\infty,
$$
which diverges at its lower limit. Uh-oh.
The most common interpretation of this divergence is to note that we never actually care about the absolute value of the potential energy; we only care about the differences between potential energies of various configurations. The quantities $U_1$ and $U_2$ don't depend on the location of the other charge; so we can view them as a fixed amount of energy that each charge carries around with it somehow. We are free to reset our "zero" for potential energy so that the potential energy of the system goes to zero as their separation $r \to \infty$, by subtracting the "constant" $U_1 + U_2$ from our definition for $U$ above; and then our new potential energy is $U_\text{int}$ by itself and everything is hunky-dory.
For most people, that resolves that, and if you're OK with the above explanation, you don't need to read the next paragraph. That said: it's still a bit odd & unsatisfying that we have these infinities running around. Really, what this calculation is telling us is the physics version of GIGO. Point charges have infinite charge density, and so we shouldn't be surprised when other important quantities (like energy) also end up being infinite when we use such ill-behaved charge distributions. If we model the charges as uniform balls with radius $R$, this whole problem never arises (though the integral $U_\text{int}$ becomes much harder to calculate exactly.) Classical electrodynamics is full of problems that arise from taking the idea of a "point charge" too seriously (I'm looking at you, Abraham-Lorentz force), and it's better to keep in the back of your mind that "point charges" are an idealization that can occasionally bite you.