The theoretical scattering cross section of electron energy loss spectroscopy (EELS) is essentially given by with the energy loss and the momentum transfer . The macroscopic dielectric function can be calculated from first principles using time-dependent density-functional theory. However, experimental EELS measurements have a finite resolution or, when operated in spatial resolution mode, yield a -integrated loss spectrum, which deviates significantly from EEL spectra calculated for specific momenta. On the other hand, integrating the theoretical spectra over is complicated by the fact that the integrand varies over several (typically six) orders of magnitude around . In this article, we present a stable technique for integrating EEL spectra over an adjustable range of momentum transfers. The important region around , where the integrand is nearly divergent, is treated partially analytically, allowing an analytic integration of the near divergence. The scheme is applied to three prototypical two-dimensional systems: monolayers of (semiconductor), hexagonal BN (insulator), and graphene (semimetal). Here, we are confronted with the added difficulty that the long-range Coulomb interaction leads to a very slow supercell (vacuum size) convergence. We address this difficulty by employing an extrapolation scheme, enabling an efficient reduction of the supercell size and thus a considerable speedup in computation time. The calculated -integrated spectra are in very favorable agreement with experimental EEL spectra.