Deriving the wave equation for electromagnetic waves I'm currently referring to the wave equation derivation given in "Introduction to Electrodynamics" by David J. Griffiths. It follows something like this: 
The electromagnetic wave equations are given by the equations:
\begin{equation}
 v^2_{ph}\nabla^2\textbf{E} = \frac{\partial^2 \textbf{E}}{\partial t^2}
\tag{1}\label{eq1}
\end{equation}
\begin{equation}
v^2_{ph}\nabla^2\textbf{B} = \frac{\partial^2 \textbf{B}}{\partial t^2}
\tag{2}\label{eq2}
\end{equation}
Where,
$v_{ph} = \frac{1}{\sqrt{\mu\epsilon}} $, is the velocity at which the wave travels in the medium 

 We have the Heaviside' form of Maxwell's equations in differential form,
$$\nabla . \textbf{E} = \frac{\rho}{\epsilon_0}$$
$$\nabla . \textbf{B} = 0$$
$$\nabla \times \textbf{E} = -\frac{\partial\textbf{B}}{\partial t}$$
$$\nabla \times \textbf{B} = \mu_0 \left(\textbf{J} + \epsilon_0\frac{\partial\textbf{E}}{\partial t}\right)$$
In a vacuum and charge-free space ($\rho = 0, \textbf{J} = 0$), these equations are:,
\begin{equation}
\nabla . \textbf{E} = 0
\tag{3}\label{eq3}
\end{equation}
\begin{equation}
\nabla . \textbf{B} = 0
\tag{4}\label{eq4}
\end{equation}
\begin{equation}
\nabla \times \textbf{E} = -\frac{\partial\textbf{B}}{\partial t}
\tag{5}\label{eq5}
\end{equation}
\begin{equation}
\nabla \times \textbf{B} = \mu_0\epsilon_0\frac{\partial\textbf{E}}{\partial t}
\tag{6}\label{eq6}
\end{equation}
We also have the vector identity for any vector function of space,
\begin{equation}
 \nabla \times \textbf{V} = \nabla.(\nabla\textbf{V}) - \nabla^2\textbf{V}
\tag{7}\label{eq7}
\end{equation}
We thus try to find the curl of $ (5) $
\begin{equation}
\nabla \times (\nabla \times \textbf{E}) = \nabla \times -\frac{\partial\textbf{B}}{\partial t} = -\frac{\partial}{\partial t}(\nabla \times \textbf{B})
\tag{8}\label{eq8}
\end{equation}
Plugging $(6)$ into $(8)$, we get
\begin{equation}
\nabla \times (\nabla \times \textbf{E}) = -\frac{\partial}{\partial t}\left(\mu_0\epsilon_0\frac{\partial\textbf{E}}{\partial t}\right) = -\mu_0\epsilon_0\frac{\partial^2\textbf{E}}{\partial t^2}
\tag{9}\label{eq9}
\end{equation}
Using the identity (7),
\begin{equation}
\nabla.(\nabla\textbf{E}) - \nabla^2\textbf{E} = -\mu_0\epsilon_0\frac{\partial^2\textbf{E}}{\partial t^2}
\tag{10}\label{eq10}
\end{equation}
Plugging (3) into (10), we get
\begin{equation}
0 - \nabla^2\textbf{E} = -\mu_0\epsilon_0\frac{\partial^2\textbf{E}}{\partial t^2}
\tag{11}\label{eq11}
\end{equation}
With a slight bit of rearrangement we can see that,
$$- \nabla^2\textbf{E} = -\mu_0\epsilon_0\frac{\partial^2\textbf{E}}{\partial t^2} $$
$$\nabla^2\textbf{E} = \mu_0\epsilon_0\frac{\partial^2\textbf{E}}{\partial t^2}$$
\begin{equation}
 \frac{1}{\mu_0\epsilon_0}\nabla^2\textbf{E} = \frac{\partial^2\textbf{E}}{\partial t^2}
\tag{12}\label{eq12}
\end{equation}
(12) and (1) correspond to the same equation under the conditions imposed. Following a similar line of reasoning for equation (6), we will get the equation;
\begin{equation}
  \frac{1}{\mu_0\epsilon_0}\nabla^2\textbf{B} = \frac{\partial^2\textbf{B}}{\partial t^2}
\tag{13}\label{eq13}
\end{equation}
which corresponds to equation $(2)$. Thus, we have obtained the electromagnetic wave equation from Heaviside' form of Maxwell's equations in differential form. 
 Is there a way to derive the same wave equations without assuming to be in a vacuum and charge-free space i.e $ \ \textbf{J} = 0, \rho = 0$?
 A: I will be going to set $\epsilon_0=\mu_0=1$. Now the Maxwell equations are: 
$$\nabla . \textbf{E} = \rho$$
$$\nabla . \textbf{B} = 0$$
$$\nabla \times \textbf{E} = -\frac{\partial\textbf{B}}{\partial t}$$
$$\nabla \times \textbf{B} =  \left(\textbf{J} + \frac{\partial\textbf{E}}{\partial t}\right)$$
 and we have the identity 
\begin{equation}
 \nabla \times (\nabla \times\textbf{V}) = \nabla(\nabla.\textbf{V}) - \nabla^2\textbf{V}
\end{equation}
Now proceeding in the same fashion as in the vacuum
$$\nabla \times (\nabla \times \textbf{E}) = \nabla \times -\frac{\partial\textbf{B}}{\partial t} = -\frac{\partial}{\partial t}(\nabla \times \textbf{B})=-\frac{\partial}{\partial t}\left(\textbf{J} + \frac{\partial\textbf{E}}{\partial t}\right)$$
while the LHS becomes:
$$\nabla(\nabla.\textbf{E}) - \nabla^2\textbf{E}=\nabla(\rho) - \nabla^2\textbf{E}$$
Rearranging RHS and LHS we get
$$\nabla^2\textbf{E}-\frac{\partial^2\textbf{E}}{\partial t^2}=\nabla\rho +\frac{\partial}{\partial t}\textbf{J}$$
In simpler terms $$\Box\textbf{E}=\textbf{C}$$
where $$\textbf{C}=\nabla\rho +\frac{\partial}{\partial t}\textbf{J}$$
Now moving to the case of $\textbf{B}$
$$\nabla \times (\nabla \times \textbf{B})=\nabla \times\left(\textbf{J} + \frac{\partial\textbf{E}}{\partial t}\right)= \nabla \times\textbf{J} + \frac{\partial}{\partial t}(\nabla\times\textbf{E})=\nabla \times\textbf{J} -\frac{\partial^2\textbf{E}}{\partial t^2}$$
as for LHS we have
$$\nabla(\nabla.\textbf{B}) - \nabla^2\textbf{B}=\nabla(0) - \nabla^2\textbf{B}$$
Rearranging RHS and LHS we get
$$\nabla^2\textbf{B}-\frac{\partial^2\textbf{B}}{\partial t^2}=-\nabla \times\textbf{J} $$
in simpler terms $$\Box \textbf{B}=\textbf{F}$$ 
where $$\textbf{F}=-\nabla \times\textbf{J}$$
So putting sources has ultimately lead to what we call inhomogeneous wave equation which is simply $$\Box f(t,\vec{x})=h(t,\vec{x})$$ same thing as in case of Laplacian and Poisson equation in chapter 3.
Bonus material(I will make an assumption of tensors): 
Maxwell equations are Lorentz covariant equation (that's how they contributed in Einstein triumph of special relativity), even when they were discovered in the era of Newtonian mechanics. Lorentz covariance another term to say a given physical quantity obey transformation law of different inertial reference frames in special relativity. 
You might have also noticed how messy it becomes taking curl and div each time in above calculation and you will see it when you compare equation of chapters 10 and 12 of Griffiths book relating $\vec{J},\rho, A_\mu$. I will provide a rough sketch of the above calculation in the light of SR.
We define a quantity called 4-vector a generalization of vectors in 4 dimensions of Minkowski space $$A_{\mu}=(V, A_x, A_y, A_z)$$
$$J_{\mu}=(\rho, J_x, J_y, J_z)$$ 
Define a quantity called electromagnetic strength tensor $$F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu$$
Maxwell equation can be recast as $$\partial^\nu F_{\mu\nu}=J_\mu$$ and $$\partial_{[\mu} F_{\nu\lambda]}=0$$
Leave the second equation aside (it's actually a tautology) let's focus on first equation expanding it in terms of $A_\mu$
$$\partial^\nu(\partial_\mu A_\nu-\partial_\nu A_\mu)=J_\mu$$
$$\partial^\nu(\partial_\mu A_\nu)-\partial^\nu(\partial_\nu A_\mu)=J_\mu$$ rearanging the terms we have 
$$\partial_\mu(\partial^\nu A_\nu)-(\partial^\nu\partial_\nu) A_\mu=J_\mu$$
Now we use the Lorentz gauge and set $\partial^\nu A_\nu$ so ultimately we're left with 
$$-(\partial^\nu\partial_\nu) A_\mu=J_\mu$$ which is nothing but $$-\Box A_\mu = J_\mu$$
which is simply the wave equation of different potentials under the presence of different sources you can recover $\vec{E}$, $\vec{B}$ from $A_\mu$. You might not have followed anything from this bonus material if it's your first encounter of,4-vectors, tensors, Einstein summation, Gauge transformation/freedom. What I actually wanted to show you was bear for the time being in the complicated mess of calculation and when you're done with chapter 12 of Griffith you'll have a different outlook at Electrodynamics as a whole.
A: The fact that the electric and magnetic fields obey wave equations of that form is a direct result of assuming that there are no charges or currents. If those assumptions are relaxed*, then in step 11 the term going to zero will in fact not be zero (note that your vector identity is written incorrectly; it should be $ \nabla \times \nabla \times \textbf{E} = \nabla(\nabla . \textbf{E}) - \nabla^2\textbf{E})$. What you end up with is:
$$
\frac{1}{\epsilon_0}\nabla \rho - \nabla^2\textbf{E} = -\mu_0 \epsilon_0 \frac{\partial^2 \textbf{E}}{\partial t^2}
$$
The presence of this extra term means that this is no longer what we consider a wave equation; in general it will not be linear, and it certainly won't have nice sinusoidal solutions. 
Also note that if you want to consider these equations not in a vacuum, if the wave is travelling in a linear, homogeneous material**, you can just replace $ \mu_0 $ and $ \epsilon_0$ with the $\mu$ and $\epsilon$ of the medium.
*For simplicity, I am assuming a time-invariant current so that it disappears in the time derivative in (8), but you could easily relax this assumption and arrive at a similar conclusion with a different final form.
** Without these assumptions, $\mu$ and $\epsilon$ will depend on space and again your equations will take on a different form with different solutions.
A: It is much easier to derive the wave equation by the Maxwell equations written in covariant form. They then read $$\partial_\mu F^{\mu\nu} = -j^\mu / \epsilon_0 ~~.$$ As $F^{\mu\nu} = \partial^mu A^\nu - \partial^\mu A^\nu$, this becomes $$\partial_\mu \partial^\mu A^\nu - \partial_\mu \partial^\nu A^\mu = -j^\mu/\epsilon_0 ~~.$$ In the Lorenz gauge , $\partial_\mu \partial^\mu A^\nu - \partial_\mu j^\mu = 0$, this becomes the wave equation for the potential $$\partial_\mu \partial^\mu A^\nu = j^\mu/\epsilon_0 ~~.$$ If you want you can find $E$ and $B$ directly from $A^\mu$.
