First of all, I think you're missing a $-\textbf{J}.\textbf{E}$ term in the RHS of your final expression. The rest of the expression looks fine.
I present here some general guidelines on how to approach this derivation. As per the homework guidelines of stackexchange I will not provide all the steps. Others are welcome to correct me on this if I have not completely understood the guidelines. I understand that solving coupled equations using vector calculus can be overwhelming and error-prone. Therefore, I will provide “anchor points,” which are nothing but validation steps that you are heading in the right direction. My TA used this technique. If you are getting some horrible terms, which I failed to mention, then it is “probably” time to step back and recheck your calculations. The reason I say “probably” is because it is possible that you come up with an alternate derivation. I think the one I worked out is the simplest one. Here it is:
You can observe that the identity is nothing but:
$$\nabla.\textbf{S}=(\nabla \times \textbf{E}).\textbf{H}-(\nabla \times \textbf{H}).\textbf{E}$$
The two terms on the RHS of the above equation should give a hint as to which equations you should manipulate first. Where in the above list can you find $\nabla \times \textbf{E}$ or $\nabla \times \textbf{H}$? You’ll need to bring those equations in that form first, i.e. the form $(\nabla \times \textbf{E}).\textbf{H}$ and $(\nabla \times \textbf{H}).\textbf{E}$ and then substitute their modified RHS in the identity. After all these manipulations, say you have obtained a form (*). You can observe that the RHS of (*) has time derivatives. You can now start seeing that (*) is closer to the form that you want. You will now need to use $\textbf{D}= \epsilon_0 \textbf{E}+\textbf{P}$ and $\textbf{B}= \mu_0 (\textbf{H}+\textbf{M})$ in (*) in the time derivatives. Yes, you are missing the latter in the above list. After a little bit of simplification, you will need to obtain a contracted form. I will show one (of the two):
$$\textbf{E}.\frac{\partial \textbf{E}}{\partial t} = \frac{\partial}{\partial t}\left(\frac{1}{2}\textbf{E}^2\right)$$
After performing a similar manipulation for the magnetic field term, you will obtain the desired expression (with the missing $-\textbf{J}.\textbf{E}$).
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).
Now electromagnetic field 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$} \;. $$ This definitely implies the conservation of energy and moreover, the equation which OP wrote in his query is the offshoot of this same continuity equation.
Energy density and Poynting Vector:
We can write the differential form of the continuity equation above as:
$$-\frac{\partial u}{\partial t}= \textrm{div}\,\mathbf S + \mathbf E\cdot \mathbf J$$ (since, all the derivation is done in Cartesian coordinates, then $\textrm{div}\equiv \mathbf \nabla$).
Or, we can write $$\mathbf E\cdot \mathbf J= -\frac{\partial u}{\partial t}-\mathbf E\cdot \mathbf J\tag 1$$
In order to find out what $u$ and $\bf S$ actually are, it is assumed that they solely depend on the fields $\bf E$ and $\bf B\;.$
Now, using $$\mathbf J = \epsilon_oc^2\,\mathbf \nabla\times\mathbf B- \epsilon_0\frac{\partial\mathbf E}{\partial t}\,,$$ we can write $\mathbf E\cdot \mathbf J$ as $$\mathbf E\cdot \mathbf J= \epsilon_oc^2\,\mathbf E\cdot (\mathbf \nabla\times\mathbf B)- \epsilon_0\mathbf E\cdot\frac{\partial\mathbf E}{\partial t}\;.$$
Now, \begin{align}\mathbf E\cdot (\mathbf \nabla\times\mathbf B)&= (\mathbf \nabla\times\mathbf B)\cdot \mathbf E\\ &= \mathbf \nabla\cdot (\mathbf B\times\mathbf E)+ \mathbf B\cdot (\mathbf \nabla \times \mathbf E)\;.\end{align}
Therefore, \begin{align}\mathbf E\cdot \mathbf J&=\epsilon_oc^2\,\mathbf \nabla\cdot (\mathbf B\times\mathbf E)+\epsilon_oc^2\,\mathbf B\cdot (\mathbf \nabla \times \mathbf E)-\epsilon_0\mathbf E\cdot\frac{\partial\mathbf E}{\partial t}\\ &=\epsilon_oc^2\,\mathbf \nabla\cdot (\mathbf B\times\mathbf E)+\epsilon_oc^2\,\mathbf B\cdot (\mathbf \nabla \times \mathbf E)-\frac{\partial}{\partial t}\,\left(\frac12 \epsilon_0\,\mathbf E\cdot \mathbf E\right)\\ &=\mathbf \nabla\cdot (\epsilon_oc^2\,\mathbf B\times\mathbf E)+ \epsilon_oc^2 \mathbf B\cdot \left(-\frac{\partial\mathbf B}{\partial t}\right)-\frac{\partial}{\partial t}\,\left(\frac12 \epsilon_0\,\mathbf E\cdot \mathbf E\right)\\ &=\mathbf \nabla\cdot (\epsilon_oc^2\,\mathbf B\times\mathbf E) - \frac{\partial}{\partial t}\left(\frac{\epsilon_oc^2}{2}\,\mathbf B\cdot \mathbf B+ \frac{\epsilon_0}2 \,\mathbf E\cdot \mathbf E\right) \tag 2\end{align}
Comparing $(1)$ and $(2)\,,$ we get $$u = \frac{\epsilon_0}2 \,\mathbf E\cdot \mathbf E+\frac{\epsilon_oc^2}{2}\,\mathbf B\cdot \mathbf B\\ \mathbf S=\epsilon_oc^2\,\mathbf E\times\mathbf B \;.$$
This vector $\bf S\,,$ that is the energy flux vector we considered at the beginning of the making up of the continuity equation, is called the Poynting Vector.
The whole derivation is based on the continuity equation, which is the mathematical expression of conservation of energy.
OP can conclude the equation he wrote in the question by converting the differential form of equation $(2)$ into integral form; using the definition of $u$ and $\bf S$ derived above and lastly using Gauss's theorem.
References:
$\bullet$ Lectures on Physics by Feynman, Leighton, Sands.
you're saying that energy conservation is needed as an independent equation, in order to deduce that the Poynting vector gives energy flux. I don't think this is true. I think the energy continuity equation can be shown as a consequence of Maxwell's equations.
No! Never, ever I've made such statement. I've answered to your query I know the books say it represents energy flux, but how do you prove it represents energy flux?; that's it. What I've done is to clearly present the continuity equation before OP and have tried to clearly interpret each term and how they together imply conservation of energy. The continuity equation can indeed be derived from Maxwell's equation and this is what I've done to define $u$ and $\bf S\;.$ I'm really pretty upset on the allegation OP made.
But this is not a proof of conservation of energy.
It's indeed a proof of conservation of energy.
Contrast this with conservation of charge within electromagnetism. Conservation of charge is a consequence of Maxwell's equations, it does not have to be assumed independently.
Okay, let's fix this thing. Conservation of charge implies this continuity equation(this can be derivable from Maxwell's equations):
$$-\frac{\partial \rho_V}{\partial t}= \mathbf \nabla \mathbf J\;.$$
This suggests conservation of energy would look like:
$$-\frac{\partial u}{\partial t}= \mathbf \nabla \mathbf S\;.$$
But that is incomplete as it is the total energy and not just field energy which is conserved; the interaction of field with matter has to be taken into consideration.
is energy conservation a consequence of Maxwell's equations or something assumed independently from it?
Yes. The continuity equation meant for conservation of energy is derived from Maxwell's equation; that's how the definition of $u$ and $\bf S$ have been computed. OP must be confused to see in my answer first the continuity equation as if I'm saying continuity equation should be considered independently in order to consider conservation of energy. No! That was intentionally done in order to conceive that the Poynting Theorem indeed implies conservation of energy.
Best Answer
It is not an answer, but more a hint. From school I remember the simple problem which forces to accept the concept of a Poynting vector. If one considers two charged particles moving in perpendicular directions and write energy/momentum conservation for this system, the solution will contain Poynting vector explicitly. In some simple cases (like particles moving in ONE direction or non-interacting particles moving in static field) you may obviously skip this concept.