one matrix-vector multiplication per time step. The time integration can be truncated due to exponential decay factors and it determines a scheme nearly independent from the dimension *N* (number of pair states) of the system. Two further advantages come from the matrix-vector multiplication scheme, i.e., an $O(N^2)$ dependence on operations count and the possibility to effectively distribute the multiplications on several processors of a parallel computer.[@smith] Applications of this method to the calculation of the $\omega$-dependent dielectric function of bulk and surface systems have been successfully performed with good comparison either with experimental results and with outcomes of the matrix diagonalization scheme. Moreover, a further speed-up of the calculations could be obtained combining this method with the use of a model dielectric function,[@smith; @cappellinissc92; @cappelliniPRB93] which, however, has not been employed in the current calculations which involve no models at all. This scheme has the big advantage that it never requires the explicit diagonalization of the Hamilton matrix, because it needs only to evaluate the Hamilton exciton matrix times the current wave function for every time-step. Massive CPU-time saving and improved scaling of the computational work with respect to system size result. On the other hand, we should admit that 