## 3 Numerical Results  
The aim of this section is the validation of the method proposed in Section 2. The approximate solution obtained by truncating a Fredholm second-kind matrix equation converges to the exact solution of the problem as the truncation order increases. In principle, the accuracy of the solution is limited only by the machine precision; in practice, it depends on the accuracy in the numerical evaluation of the coefficient matrix. It will be shown in the following that the solution proposed in this paper, together with the analytical technique shown in [40] devised to efficiently evaluate the 1-D improper integrals of oscillating and algebraic decaying functions filling the coefficient matrix, makes it possible to accurately reconstruct the solution with very low CPU time and memory occupation. It is worth noting that all the simulations are performed by means of an in-house software code running on a laptop equipped with an Intel Core i7-10510U 1.8 GHz, 16 GB RAM. The rate of convergence will be analyzed by introducing the following relative computation error as a function of the truncation order of the matrix equation:  
\[\text{err}_{N}(M)=\sqrt{\sum_{n=-N+1}^{N-1}\left\|\mathbf{x}_{M+1}^{(n)}- \mathbf{x}_{M}^{(n)}\right\|^{2}/\sum_{n=-N+1}^{N-1}\left\|\mathbf{x}_{M}^{(n )}\right\|^{2}}, \tag{14}\]  
where \(2N-1\) is the number of cylindrical harmonics considered according to the estimation formula in [51], \(\|\cdot\|\) denotes the Euclidean norm, and \(\mathbf{x}_{M}^{(n)}\) is the vector related to the first \(M\) expansion coefficients of the general harmonics of the unknown functions. Henceforth, it will be supposed that \(z_{1}=0\) and \(z_{2}=-d\) (so that the distance between the disks could be changed by simply moving disk 2 along the \(z\) axis), \(R_{1}=100\)\(\Omega\) and \(R_{2}=200\)\(\Omega\), \(|\)\(E_{0}|=1\) V/m, \(\theta_{0}=30^{\circ}\), \(\phi_{0}=0^{\circ}\), and transverse electric (TE) incidence with respect to the z axis. In Figure 2, \(\text{err}_{N}(M)\) is plotted for different values of the radius (\(a=\lambda/2,\lambda,2\lambda\)) and the distance between the disks (\(d/a=0.1,1,10,\infty\)). As can be clearly seen, the error is substantially independent of the distance between the disks and it depends slightly on the radius of the disks in all the examined cases. It is apparent as the convergence is really very fast. Indeed, in the worst case examined (\(a=2\lambda\) and \(d/a=0.1\)), a relative error less than \(1\%\) is obtained for \(M=6\) and \(N=15\) with a computation time (\(t\)) of only 9 s. Moreover, more accurate solutions can be quickly obtained. Indeed, a relative error less than \(0.1\%\) and \(0.01\%\) is obtained for at most \(M=11\), \(t=25\) s and \(M=22\), \(t=65\) s, respectively. In the examples examined in the following \(M\) will be chosen in order to guarantee a relative error less than \(1\%\). In Figure 3, for the sake of completeness, the reconstructed behavior of the components of the effective surface current densities for \(d/a=1\) and for different values of the radius (\(a=\lambda/2,\lambda,2\lambda\)) is shown. Moreover, in Figure 4, the behavior of the components of the current on the disk 1, plotted for \(a=\lambda\) and for different values of the distance between the disks (\(d/a=0.1,1,10,\infty\)), clearly shows that the proposed method tends to reconstruct the solution obtained in the case of a single disk (\(d/a=\infty\)) as the distance between the disks increases. Figure 2: Relative computation error for a two stacked thin resistive disks resonator with different values of the radius of the disks (\(a=\lambda/2,\lambda,2\lambda\)) and of the distance between the disks (\(d/a=0.1,1,10,\infty\)). \(R_{1}=100\)\(\Omega\), \(R_{2}=200\)\(\Omega\), \(|\)\(\underline{E}_{0}|=1\) V/m, \(\theta_{0}=30^{\circ}\), \(\phi_{0}=0^{\circ}\), and TE incidence: (**a**) \(d/a=0.1\); (**b**) \(d/a=1\); (**c**) \(d/a=10\); (**d**) \(d/a=\infty\).