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}$).