A numerical scheme for solving the nonlinear Heisenberg-Euler equations in up to three spatial dimensions plus time is presented and its properties are discussed. The algorithm is tested in one spatial dimension against a set of already known analyti
cal results of vacuum effects such as birefringence and harmonic generation. Its power to go beyond analytically solvable scenarios is demonstrated in 2D. First parallelization scaling tests are conducted in 3D.
