the scattered and/or transmitted fields. Another, more realistic, approach is to use a current source ${\mathbf{J}}(t)={\mathbf{J}}(\br,t)$. Instead of Eq. (\[TDME\]) we have $${\frac{\partial}{\partial t}}\Psi(t)={\mathcal H}\Psi(t) - {\mathcal J}(t),
\label{TDMEJ}$$ where $\Psi(t)=(\bX(t),\bY(t))^T$ and ${\mathcal J}(t)=(0,{\mathbf{J}}(t))^T$ represents the source term. The formal solution of Eq. (\[TDMEJ\]) is given by $$\Psi(t) = e^{t{\mathcal H}}\Psi(0)-\int_0^t
e^{(t-u){\mathcal H}}{\mathcal J}(u)\,du,
\label{formalJ}$$ showing that we can simply re-use one of the unconditionally stable algorithms to compute the second term in Eq. (\[formalJ\]). In practice, for a time-step $\tau$, we update $\Psi(t)$ according to $$\Psi(t+\tau) = e^{\tau{\mathcal H}}\Psi(t)-\int_t^{t+\tau}
e^{(t+\tau-u){\mathcal H}}{\mathcal J}(u)\,du.
\label{Psitau}$$ A standard quadrature formula can be used to compute the integral over $u$ [@Michielsen98]. When a current source is present we take as the initial condition $\Psi(0)=0$.

Time-domain algorithms can also be used to compute the eigenvalues of $H$, the discretized form of ${\mathcal H}$. In general $H$ is a (very) large matrix, usually too large to be stored in memory. If only a few, well-separated eigenvalues of $H$ are required sparse-matrix techniques can be used to compute these eigenvalues [@GOLUB; @WILKINSON]. However, if one is interested in global features of the distribution of eigenvalues, i.e. if we want to determine [*all*]{} eigenvalues, time-domain algorithms offer 