## 3 Results
### System Tests  
We perform here two kinds of tests: accuracy tests (to debug possible losses of accuracy within the code) and error convergence tests. As a unit cell, we use a vacuum cube with dimensions \([0,1]\times[0,1]\times[0,1]\) m. Regarding finite elements, we have used tetrahedra with second-order basis functions from [33]. We employ GiD [32] to generate the mesh and show the results of the simulations and employ Fortran as the programming language (following an object-oriented paradigm). The mesh that we use here is composed of 8324 tetrahedra for the unit cell, which produces a problem of 53,124 unknowns for the unit cell. We have considered a grid of \(3\times 3\) unit cells. For the accuracy test, we want to detect if all the polynomials that belong to the space of functions we have used to numerically approximate the solution yields a numerically zero output. We test each of the monomials (three-dimensional polynomial vectors) whose linear combination yields the basis functions \(N_{i}\) of (11). Specifically, the monomials involved in these basis functions are the following (please refer to [33]):  
\[ \begin{split}\left\{\begin{bmatrix}1\\ 0\end{bmatrix},\begin{bmatrix}x\\ 0\end{bmatrix},\begin{bmatrix}y\\ 0\end{bmatrix},\begin{bmatrix}z\\ 0\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}0\\ 1\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}0\\ 2\\ 0\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},\dots\\ \ldots,\begin{bmatrix}y^{2}\\ -xy\\ -y\end{bmatrix},\begin{bmatrix}0\\ -yz\\ 2\end{bmatrix},\begin{bmatrix}-xy\\ 2\end{bmatrix},\begin{bmatrix}-xz\\ 0\\ 0\end{bmatrix},\begin{bmatrix}2\\ 0\\ 2^{2}\end{bmatrix},\begin{bmatrix}0\\ -xz\\ 0\end{bmatrix},\begin{bmatrix}0\\ -xz\\ 0\end{bmatrix},\begin{bmatrix}yz\\ -xz\\ 0\end{bmatrix}\right\}\end{split}  \tag{20}\]  
We show in Table 1 the value of \(\mathbf{e}_{\text{field}}\) from (15a) for all the monomials within the space of basis functions for Test 4 from Section 2.3.2. However, we have followed the same TDD approach suggested in Section 2.3.2: i.e., we have obtained good results with Test 1, then we use Test 2 until we obtain good results, and so on. For the sake of brevity, we do not show the value for all the polynomials but rather the maximum value in Table 2. Moreover, for brevity reasons, we only show the values of \(\mathbf{e}_{\text{field}}\), although similar results are obtained with \(\mathbf{e}_{\text{rot}}\). We obtain numerically equivalent values for the three algorithms since we follow a mathematical approach, and we should obtain the same solution for all of them. Furthermore, we use a double-precision accuracy in our code so the minimum resolution is \(10^{-16}\). The three approaches yield the same numerical solution, whereas the differences are due to the different operations applied to obtain the final solution. The accuracy of the results obtained (in the order of \(10^{-12}\)) is due to the numerical noise, and it sets the numerical resolution for this problem. Second, we perform a test on error convergence, specifically using a plane wave as the solution. We use a plane wave because it is an electromagnetically relevant example of a smooth function for which theoretical results concerning the convergence of the error exist so we can compare with it. Specifically, for smooth functions and the second-order finite element basis functions of [33], the theory predicts the error to decrease as \(O(h^{2})\)[34], i.e., the slope in a log-log plot of the error versus the discretization size \(h\) should be ideally equal to 2. The expression of the electric field of a linearly polarized (in \(\hat{\mathbf{\theta}}+\hat{\mathbf{\phi}}\)) plane wave incident from \((\theta,\phi)\) direction is:  
\[\mathbf{E}=E_{0} \begin{cases}\cos(\theta)\cos(\phi)-\sin(\phi)\\ \cos(\theta)\cos(\phi)+\cos(\theta)\\ -\sin(\theta)\end{cases}