  We present a certifiable algorithm to calculate the eigenvalue density
function -- the number of eigenvalues within an infinitesimal interval -- for
an arbitrary 1D interacting quantum spin system. Our method provides an
arbitrarily accurate numerical representation for the smeared eigenvalue
density function, which is the convolution of the eigenvalue density function
with a gaussian of prespecified width. In addition, with our algorithm it is
possible to investigate the density of states near the ground state. This can
be used to numerically determine the size of the ground-state energy gap for
the system to within a prespecified confidence interval. Our method exploits a
finitely correlated state/matrix product state representation of the propagator
and applies equally to disordered and critical interacting 1D quantum spin
systems. We illustrate our method by calculating an approximation to the
eigenvalue density function for a random antiferromagnetic Heisenberg model.
