We present a method to considerably improve the numerical performance for
solving Eliashberg-type coupled equations on the imaginary axis. Instead of the
standard practice of introducing a hard numerical cutoff for treating the
infinite summations involved, our scheme allows for the efficient calculation
of such sums extended formally up to infinity. The method is first benchmarked
with isotropic Migdal-Eliashberg theory calculations and subsequently applied
to the solution of the full-bandwidth, multiband and anisotropic equations
focusing on the FeSe/SrTiO3 interface as a case study. Compared to the
standard procedure, we reach similarly well converged results with less than
one fifth of the number of frequencies for the anisotropic case, while for the
isotropic set of equations we spare approximately ninety percent of the
complexity. Since our proposed approximations are very general, our numerical
scheme opens the possibility of studying the superconducting properties of a
wide range of materials at ultra-low temperatures