## 6 Numerical Methods  
This section is devoted to brief descriptions of the numerical methods used to approximate the system (14)-(19). Other techniques for similar problems were employed in [17; 22]. Preliminary computations based on this specific cylindrical setting were tried in [23]. Our approach sticks to the framework of low-order finite difference schemes. More accurate numerical methods may be applied in future developments. For the reasons we are going to describe, we soon discarded the idea of applying the popular [PERSON] scheme (see, e.g., [24; 25]). Consider the two Maxwell equations \(\partial\mathbf{E}/\partial t=c^{2}\
abla\times\mathbf{B}\) and \(\partial\mathbf{B}/\partial t=-\
abla\times\mathbf{E}\) relative to a region of the void space in absence of sources of any kind. Let us place ourselves in a Cartesian reference frame \((x,y,z)\). [PERSON]'s algorithm does not explicitly implement the equations \(\
abla\cdot\mathbf{E}=0\) and \(\
abla\cdot\mathbf{B}=0\). Indeed, the space grids are structured in such a way to reproduce only the action of the curl operator. As a consequence of this choice, the divergence-free conditions are implicitly satisfied. In a simplified context, we may assume that the fields are of the form \(\mathbf{E}=(E_{x},E_{y},0)\), \(\mathbf{B}=(0,0,B_{z})\) and do not depend on the variable \(z\). The scheme uses a system of two staggered grids of step \(h\). The fields are approximated on such grids, in order to ensure that every \(\mathbf{E}\) component is surrounded by four \(\mathbf{B}\) components, and vice versa. Time discretization is performed by fractional steps, by alternating the evaluation of \(\mathbf{E}\) and \(\mathbf{B}\). In Cartesian coordinates with grid points \((x_{i},y_{j})\), the procedure reads as follows:  
\[(E_{x})_{i-\frac{1}{2},j+\frac{1}{2}}^{k+\frac{1}{2}}=(E_{x})_{i-\frac{1}{2}, j+\frac{1}{2}}^{k-1}+\frac{c^{2}\Delta t}{h}\Big{(}(E_{x})_{i-\frac{1}{2},j+1}^{k}-(E_ {x})_{i-\frac{1}{2},j}^{k}-(E_{y})_{i,j+\frac{1}{2}}^{k}+(E_{y})_{i-1,j+\frac{ 1}{2}}^{k}\Big{)}\]  
\[(E_{x})_{i-\frac{1}{2},j}^{k}=(E_{x})_{i-\frac{1}{2},j}^{k-1}+\frac{c^{2} \Delta t}{h}\Big{(}(B_{z})_{i-\frac{1}{2},j+\frac{1}{2}}^{k-\frac{1}{2}}-(B_{ z})_{i-\frac{1}{2},j-\frac{1}{2}}^{k-\frac{1}{2}}\Big{)}\]\[(E_{y})_{i,j+\frac{1}{2}}^{k}=(E_{y})_{i,j+\frac{1}{2}}^{k-1}-\frac{c^{2}\Delta t}{ h}\Big{(}(B_{2})_{i+\frac{1}{2},j+\frac{1}{2}}^{k-\frac{1}{2}}-(B_{2})_{i-\frac{1}{2}, j+\frac{1}{2}}^{k-\frac{1}{2}}\Big{)} \tag{26}\]  
where the index \(k\) refers to time. The numerical divergence of the electric field is calculated as follows:  
\[\rho_{i,j}^{k}=\frac{1}{h}\Big{(}(E_{x})_{i+\frac{1}{2},j}^{k}-(E_{x})_{i-\frac {1}{2},j}^{k}+(E_{y})_{i,j+\frac{1}{2}}^{k}-(E_{y})_{i,j-\frac{1}{2}}^{k}\Big{)}. \tag{27}\]  
It is easy to check that \(\rho_{i,j}^{k}=\rho_{i,j}^{k-1}\), for any \(k\geq 1\), which means that such a quantity is preserved during the time iterations. Many approximation schemes for the electrodynamics equations trust the fact that the divergences will not be modified during the evolution. Unfortunately, this may be true only for isolated regions of space. When any wave interacts with other objects, for example, incompatibilities grow by putting together the evolution equations, the divergence-free conditions and the boundary constraints. This fact has been made evident in [26; 27], where it is shown that these restrictions are too many, so that they cannot hold together at the same time. By allowing \(\
abla\cdot\mathbf{E}\) to be different from zero we can come out from this dead end, thereby justifying the adoption of the consistent model we are examining in this paper. At this point, it should be clear that [PERSON]'s scheme and similar ones are not appropriate for our study. We then adopt a classical explicit (first-order) Euler method for advancing in time. The space derivatives in (14)-(16) are approximated by centered (second-order) differences.