Although this is a standard derivation, you frequently don't see it in introductory electromagnetism courses, maybe because those courses shy away from the heavy use of vector calculus. Here's the usual approach. We'll find a wave equation from Maxwell's equations.
Start with
$\nabla \times \vec{E} = -\frac{\partial\vec{B}}{\partial t}$.
Take a partial derivative of both sides with respect to time. The curl operator has no partial with respect to time, so this becomes
$\nabla \times \frac{\partial\vec{E}}{\partial t} = -\frac{\partial^2\vec{B}}{\partial t^2}$.
There's another of Maxwell's equations that tells us about $\partial\vec{E}/\partial t$.
$\nabla \times \vec{B} = \mu_0\epsilon_0\frac{\partial \vec{E}}{\partial t}$
Solve this for $\partial\vec{E}/\partial t$ and plug into the previous expression to get
$\nabla \times \frac{(\nabla \times \vec{B})}{\mu_0\epsilon_0} = -\frac{\partial^2 \vec{B}}{\partial t^2}$
the curl of curl identity lets us rewrite this as
$\frac{1}{\mu_0 \epsilon_0}\left(\nabla(\nabla \cdot \vec{B}) - \nabla^2\vec{B}\right) = -\frac{\partial^2 \vec{B}}{\partial t^2}$
But the divergence of the magnetic field is zero, so kill that term, and rearrange to
$\frac{-1}{\mu_0 \epsilon_0}\nabla^2\vec{B} + \frac{\partial^2 \vec{B}}{\partial t^2} = 0$
This is the wave equation we're seeking. One solution is
$\vec{B} = B_0 e^{i (\vec{x}\cdot\vec{k} - \omega t) }$.
This represents a plane wave traveling in the direction of the vector $\vec{k}$ with frequency $\omega$ and phase velocity $v = \omega/|\vec{k}|$. In order to be a solution, this equation needs to have
$\frac{\omega^2}{k^2} = \frac{1}{\mu_0\epsilon_0}$.
Or, setting $v = 1/\sqrt{\mu_0\epsilon_0}$
$\frac{\omega}{k} = v$
This is called the dispersion relation. The speed that electromagnetic signals travel is given by the group velocity
$\frac{d\omega}{d k} = v$
So electromagnetic signals in a vacuum travel at speed $c = 1/\sqrt{\mu_0\epsilon_0}$.
Edit
You can follow the same steps to derive the wave equation for $\vec{E}$, but you will have to assume you're in free space, i.e. $\rho = 0$.
Edit
The curl of the curl identity was wrong, there's a negative number in there
Best Answer
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}$).