Intuition for energy density of electromagnetic fields I have the following equations:
$U = \frac{1}{2}(\mathbf{E}\cdot\mathbf{D} + \mathbf{B}\cdot\mathbf{H})$ and
$\mathbf{S} = \mathbf{E} \times \mathbf{H}$
This is from the equation:
$$\frac{\partial U}{\partial t} + \mathbf{\nabla} \cdot \mathbf{S} = - \mathbf{J} \cdot \mathbf{E}$$
Although it is easy to see that $U$ has the units of energy per unit volume, I would like some intuition behind WHY it represents the energy density of electromagnetic fields.
 A: Take a spherical shell of charge with radius $R$ and total charge $Q$. How much work is needed to assemble this collection? Now shrink it to radius $r < R$. How much energy is needed to make it now?
Next, ask yourself, where is that energy stored? Although when we teach about potential energy we ascribe it as a property of the particle, that can't be where the energy is stored.
Consider the following experiment. Start with a positive charge $q$ at the origin. Bring another positive charge $Q$ in from infinity to the point $(x,0,0)$. How much energy is stored in the system? Obviously, $qQ/(4\pi\epsilon_0 |x|)$. Given how I described how to construct this system, you might think that the energy is stored in the charge $Q$, but if we move $q$ to infinity you have to ask the question: when and how did the energy leave $Q$? Keep in mind: the energy has to move at the speed of light or less to satisfy special relativity.
The answer is this: it was never in $Q$ to begin with. The energy was stored in the change to the electric field that comes from the interaction of the two point charges. Specifically,
\begin{align}
  \mathbf{E}_{\text{tot}} &= \mathbf{E}_1 + \mathbf{E}_2
\end{align}
The energy stored in the electric field of a point charge is infinite, but it is also unchanging, so we have to only consider how the energy density is changed by the interaction
\begin{align}
  \Delta U_{1,2} &= \int \mathrm{d}^3x \epsilon_0 \mathbf{E}_1\cdot \mathbf{E}_2.
\end{align}
You can construct similar arguments based on electric currents and/or particle spins for the magnetic field, but they're a bit more technically involved.
Fun fact: the potential energy density of the electromagnetic field is $\frac{1}{2}\mathbf{B}\cdot\mathbf{H}$ and the kinetic energy density $\frac{1}{2}\mathbf{E}\cdot\mathbf{D}$.
A: You got that equation wrong, it should be
$$
- \mathbf{J} \cdot \mathbf{E} = \frac{\partial U}{\partial t} + \nabla \cdot \mathbf{S}
$$
where
$$
U = \frac{1}{2}\epsilon_0 E^2 +\frac{1}{2\mu_0}B^2,
$$
$$
\mathbf S = \mathbf E\times \mathbf B/\mu_0~~~(*).
$$
This is so-called Poynting's theorem (differential version). It follows directly from Maxwell's equations alone, provided the fields and sources are not too singular. So by itself this does not imply anything about mechanical concepts such as work or energy.
However, this equation contains the term
$$
- \mathbf{J} \cdot \mathbf{E}
$$
which often can be interpreted (e.g. current in a wire) as minus rate of work of electric force on mobile charges (in unit volume). Integrating this equation over some fixed space region $V$, we get integral version of the Poynting theorem, which is
$$
\int_V - \mathbf{J} \cdot \mathbf{E} ~dV= \int_V \frac{\partial U}{\partial t}dV + \int_V \nabla \cdot \mathbf{S} ~dV.
$$
Using the Gauss theorem, we can transform the last volume integral into surface integral and write
$$
\int_V - \mathbf{J} \cdot \mathbf{E}~dV = \int_V \frac{\partial U}{\partial t}dV + \oint_\Sigma d\boldsymbol{\Sigma}\cdot \mathbf{S}.
$$
If the left-hand side can be interpreted as minus work done by electric force, then we can interpret the right-hand side as
$$
\text{rate of change of EM energy inside V} + \text{rate of loss of EM energy from V through the boundary surface $\Sigma$}.
$$
From this, it is natural and the simplest possibility to define density of EM energy to be given by $U$, and energy flux density by $\mathbf S$, as given in (*). That is, however, not the only possibility for density of energy and density of energy flux, because any pair of quantities $U',\mathbf S'$:
$$
U' = U+f,
$$
$$
\mathbf S' = \mathbf S + \mathbf K,
$$
where $f, \mathbf K$ are any field functions that obey the equation
$$
\partial_t f + \nabla \cdot \mathbf K = 0,
$$
would be just as valid for density of energy and density of energy flux. But for simplicity, most commonly used expressions are those in (*) (the Poynting expressions).
Poynting's equation fails at point particles, so in case the particles are points, one has to define EM energy and momentum via different expressions, e.g. like Frenkel did [1], for a short explanation see [2].
[1] J. Frenkel, Zur Elektrodynamik punktfoermiger Elektronen, Zeits. f. Phys., 32, (1925), p. 518-534. https://doi.org/10.1007/BF01331692
[2] R. C. Stabler, A Possible Modification of Classical Electrodynamics, Physics Letters, 8, 3, (1964), p. 185-187. https://doi.org/10.1016/S0031-9163(64)91989-4
A: Look up Griffiths introduction to electrodynamics page 92, for the derivation of electrostatic potential energy which is the work done  against the electric field to assemble a charge distribution measured from some reference location( taken mostly as infinity)
Look up page 317 for the magentic field energy density derivation in magnetostatics, which represents the work done to accelerate charges against the induced electric field created due to the changing magnetic field when charges are accelerating. ( energy is stored in the B field and is released when the B field dissappears via induction)
these are derived in ELECTROSTATICS and magnetoSTATICS however they are also derived more generally in poyntings theorem.
Poyntings theorem can be derived without any preconceived knowledge about energy conservation on page : 346
This involves relating the work done on a general charge distribution  due to the lorentz force, and then substituting maxwells equations into this .
It's the derivation in electro and magnetostatics that gives you the insight , but derivation of poyntings theorem gives it more generally
