Field is not any abstract fairy-tale 'mathematical fiction' concept; it has momentum and energy density; it exchanges momentum and energy with charged entities.
Energy conservation is a local process which evidently implies electromagnetic field between two interacting charges must mediate the energy and momentum exchange between the charges and hence must have energy density and momentum.
Let the electromagnetic field has $u$ as its energy density (amount of energy per unit volume in the field) and let $\bf S$ represents the energy flux- the amount of energy per unit time flowing across a unit area perpendicular to the flow).
Electromagnetic field, as said at the outset, can interact with matter and do work on them; so this energy interaction must be considered when discussing energy conservation.
The field energy inside a volume $V$ is $\displaystyle \int_V u\,\mathrm dV\;.$ Amount of energy flowing out of the volume $V$ is given by $\displaystyle \int_\Sigma \mathbf S\cdot n\,\mathrm da \;.$
Now, the work done per unit time by the field on the matter inside the volume $V$ is given by $\displaystyle \int _V Nq(\mathbf E + \mathbf v\times \mathbf B)\cdot \mathbf v\,\mathrm dV$ where $N$ is the number of particles per unit volume; this can be written as \begin{align}\int _V Nq(\mathbf E + v\times \mathbf B)\,\mathrm dV& = \int_V Nq\mathbf E\cdot v\,\mathrm dV\\ &= \int_V \mathbf E\cdot \underbrace{(Nq\mathbf v)}_\textrm{current-density}\,\mathrm dV\\ &= \int_V E\cdot \mathbf J\,\mathrm dV \end{align}
So, the continuity equation is written thus: $$\underbrace{-\frac{\partial}{\partial t}\int_V u\,\mathrm dV}_\textrm{rate of change of energy inside volume $V$}=\underbrace{\int_\Sigma \mathbf S\cdot \mathbf n\,\mathrm da}_\textrm{amount of $\mathbf{field\, energy}$ flowing out of volume $V$ per unit time} + \underbrace{\int_V \mathbf E\cdot \mathbf J\,\mathrm dV}_\textrm{work done per unit time by the field on the matter inside volume $V$} \;. $$
After some algebra, it can be inferred that
$$u = \frac{\epsilon_0}2 \,\mathbf E\cdot \mathbf E+\frac{\epsilon_oc^2}{2}\,\mathbf B\cdot \mathbf B\;.$$
It boils down to much simpler interpretation in the case of statics as aired by Jahan Claes.
The energy $U$ is the potential energy associated with configuration of the charges of the system.
If you keep three charges very, very far from each other, then the potential energy of the system is very effectively zero.
If you keep three charges very, very far from each other, then the potential energy of the system is very effectively zero.
But when you bring them close together to a specified coordinate, then the potential energy of the system increases from $0$ to a certain positive value $U\;.$
For a discrete system of $N$ charges, the potential energy associated with their configuration is given by
\begin{align}U&= \frac12\,\sum_{j=1}^N \, q_j\sum_{k\ne j}\,\frac1{4\pi\epsilon_0}\,\cdot \frac{q_k}{r_{jk}}\\ &= \frac12\,\sum_{j=1}^N \, q_j\,\varphi(\mathbf r_j) \end{align}
For a system of continuous charges, $$U= \frac12\int \rho(\mathbf r)\,\varphi(\mathbf r)\,\mathrm d^3\mathbf r \;.$$
After some manipulation, it can again be inferred that
$$U= \frac{\epsilon_0}2 \int\,\mathbf E\cdot \mathbf E\,\mathrm d^3\mathbf r\;.$$
This is the energy required to assemble the system of charges in that specific configuration.