The physics of waves control most of the world, in multiple forms, such as electromagnetic waves. Mathematicians and physicists have developed equations which describe the patterns in which waves evolve over time, while moving through space. Due to their partial differential form, solutions to these equations must be approximated. This study introduces a new numerical scheme to perform the approximation which is highly stable and computationally efficient. This numerical scheme is formulated with respect to Maxwell’s equations, employing spatial and temporal staggering to implement a fourth-order phase accuracy. It is then compared to the traditional Yee scheme and the Runge-Kutta 3 scheme in one-dimensional applications, revealing a similar accuracy to the Runge-Kutta 3 scheme while requiring less computations per time step. Simulations are then performed in the two-dimensional case. First, no boundary conditions are implemented, causing reflection at the edge of the spatial domain. Next, the simulation is conducted while employing absorbing boundary conditions, simulating wave propagation over an infinite spatial domain. These results are compared to the results of a large domain simulation, in which the wave propagation does not reach the boundaries. Comparing the simulations, it is concluded that the numerical scheme is stable and highly accurate when employing absorbing boundary conditions. Finally, the scheme is tested in two dimensions with wave propagation through nonlinear media, as opposed to the prior simulations which were performed as if in a vacuum. After performing spectral analysis on the resulting waves after a long-time domain simulation, the resulting angular frequencies match those expected from theory. Therefore, the scheme is concluded to be powerful in one-dimensional, two-dimensional, and nonlinear simulations, all while being computationally efficient.