We present a novel method for the calculation of the energy density of states D(E) for systems described by classical statistical mechanics. The method builds on an extension of a recently proposed strategy that allows the free energy profile of a canonical system to be recovered within a pre-assigned accuracy,[A. Laio and M. Parrinello, PNAS 2002]. The method allows a good control over the error on the recovered system entropy. This fact is exploited to obtain D(E) more efficiently by combining measurements at different temperatures. The accuracy and efficiency of the method are tested for the two-dimensional Ising model (up to size 50x50) by comparison with both exact results and previous studies. This method is a general one and should be applicable to more realistic model systems.