Why is FDTD derived directly from Maxwell's equations instead of the wave equation I've been wondering why the Finite Difference Time Domain Method is derived directly from Maxwell's equations and not directly from the electromagnetic wave equation (that in theory is also derived from Maxwell's equations).
So instead of doing a Yee grid for these two differential equations
\begin{gather}
\mathbf{\nabla} \times \mathbf{E} = \  - \frac{\partial\mathbf{B}}{\partial t},   \label{ec4}
\\
\mathbf{\nabla} \times \mathbf{B =}\mu\left( \mathbf{J} + \varepsilon\frac{\partial\mathbf{E}}{\partial t} \right)   ,
\end{gather}
why don't we just do a finite-difference method for this equation:
\begin{equation}
\nabla^{2}\mathbf{E} - \mu\varepsilon\frac{\partial^{2}\mathbf{E}}{{\partial t}^{2}} = \ \nabla\left( \frac{\rho}{\varepsilon} \right) + \mu\frac{\partial\mathbf{J}}{\partial t}   
\end{equation}
I know that the first two only involve first-order derivatives and the third one is second-order (because of the Laplacian) but the third one has implicit the magnetic field, so there will be no need to reorganize both the electric and magnetic field in a convenient way as Yee proposed, because we will only need to numerically approximate the electric field, so it will be possible to use the collocated grid. 
I'm 100% sure that this is not new at all, and that I am not taking into account a lot of important and basic theories (I have to admit that I'm not an expert in this subject).  I've been trying to find an answer but google hasn't been helpful so far.
 A: It is possible to solve for the electric field $\mathbf{E}$ given the wave equation you've written down.  As I hinted in my comment, there is no standard way to directly integrate equations that are second-order in time.  However, there is a standard method to deal with such equations: introduce an auxiliary variable that represents $\partial\mathbf{E} / \partial t$.  Let's call it $\mathbf{F}$.  Then, the system of equations you actually evolve is
\begin{gather}
  \frac{\partial \mathbf{E}}{\partial t} = \mathbf{F},
  \\
   \frac{\partial\mathbf{F}}{{\partial t}} = - \frac{1}{\mu\varepsilon} \left( \nabla\left( \frac{\rho}{\varepsilon} \right) + \mu\frac{\partial\mathbf{J}}{\partial t} - \nabla^{2}\mathbf{E}\right).
\end{gather}
Now, you and I know that these two equations really represent the single wave equation above, but the numerical ODE integrator that you have doesn't.  Moreover, we can actually apply theorems about the well-posedness of this system of equations, which tell us that solutions exist (and maybe that our ODE integrator is stable in certain circumstances).
But now, by decomposing your wave equation in this way, we've introduced another couple headaches to the actual problem: setting initial values and boundary conditions for $\mathbf{F}$.  On top of that, what if you want to actually understand the physics in this situation?  You need to compute $\mathbf{B}$ somehow.  And how do you do that?  Well, plugging in the known value $\mathbf{E}$, and then evolving your first Maxwell's equation.  What has this approach bought you?
