A promising time domain electromagnetics numerical method for treating the highly nonlinear problem of charge transport in electronic devices called Delaunay-Voronoi surface integration is presented. This method couples the rotational electric and magnetic fields governed by Ampere's and Faraday's laws with the electrostatic potential dictated by Poisson's equation in a simultaneous solution. Discretization of the governing equations using dual meshes and the relevant boundary conditions are presented. The engineering application details specific to electronic device simulation are treated, and an example calculation is shown to compare with an analytical solution for propagation in a waveguide. Benchmark results are presented for the rotational equations, Poisson's equation, and the complete set of electromagnetic equations. the computational domain, and electric and magnetic fields are computed through scalar and vector potentials. An IE is most often derived from Maxwell's equations using the Green's function approach [2]. The crux of the approach is deriving a Green's function that satisfies the linear Helmholtz equation subject to point sources. IEs are then built upon the principle of linear superposition. Although this method is a staple of antenna design and other radiation (free-space) applications, the Green's function method cannot be applied when the point source is inherently nonlinear and is very limited when the dielectric medium is inhomogeneous. Other methods for deriving IEs from Maxwell's equations may produce viable options for coupling to electronic transport, but they are unknown to the authors. For a more detailed discussion of IEs, please see [2] and the references therein. Time-domain differential equation methods are much better suited for representing and solving full-wave EM phenomena in electronic devices. Unlike IE methods, DE methods seek EM solutions directly from Maxwell's equations in differential form. The two most common TD-DE categories can be classified as unstructured and structured grid methods. Even though these two methods are both DE solvers, they differ greatly in their computational approach and sophistication. Unstructured grid (mesh) methods first decompose the computational domain into a collection of elements that are typically tetrahedra. The governing equations are then written in a form suitable for discretization. This amounts to writing Maxwell's equations in curl-curl form for the well-known finite element method TD [2] or adding flux terms to Maxwell's equations for the discontinuous Galerkin TD (DGTD) method [3]. Both of these methods seek a variational solution by approximating the quantities of interest with collections of basis functions, which are carefully tailored to satisfy the governing equations and boundary conditions (BCs), for example, vector finite elements. When coupling these discretization methods to charge transport, it is not clear to the authors how to construct a basis set that approximates an equation with exponentially nonlinear depend...