The mass returned to the ambient medium by aging stellar populations over cosmological times sums up to a significant fraction (20% - 30% or more) of their initial mass. This continuous mass injection plays a fundamental role in phenomena such as galaxy formation and evolution, fueling of supermassive black holes in galaxies and the consequent (negative and positive) feedback phenomena, and the origin of multiple stellar populations in globular clusters. In numerical simulations the calculation of the mass return can be time consuming, since it requires at each time step the evaluation of a convolution integral over the whole star formation history, so the computational time increases quadratically with the number of time-steps. The situation can be especially critical in hydrodynamical simulations, where different grid points are characterized by different star formation histories, and the gas cooling and heating times are shorter by orders of magnitude than the characteristic stellar lifetimes. In this paper we present a fast and accurate method to compute the mass return from stellar populations undergoing arbitrarily complicated star formation histories. At each time-step the mass return is calculated from its value at the previous time, and the star formation rate over the last time-step only. Therefore in the new scheme there is no need to store the whole star formation history, and the computational time increases linearly with the number of time-steps. ; Comment: 8 pages, 4 figures, submitted, comments welcome