We present the generalization of the recursion method of Haydock and co-workers to systems of many interacting particles. This new method has close similarities to the memory function or Mori formalism, but with some important differences. Heisenberg's equation for the time evolution of a microscopic operator is recursively transformed into a tridiagonal matrix equation. This equation resolves the operator into components corresponding to transitions of different energies. The projected spectrum of transitions has a continued fraction expansion given by the elements of the tridiagonal matrix. We show that for an appropriate choice of inner product this density of transitions obeys a generalization of the black body theorem of electromagnetism, in that it is exponentially insensitive to distant parts of the system. This implies that the projected density of transitions is computationally stable and can be calculated even in macroscopic many-body systems. We argue that the physical content of the density of transitions is determined by the nature of its singular points, such as discrete transitions, continuous spectrum, band edges and van Hove singularities.