We propose a method based on the Wang-Landau algorithm to numerically generate the spectral densities of random matrix ensembles. The method employs Dyson's log-gas formalism for random matrix eigenvalues and also enables one to simultaneously investigate the thermodynamic properties. This approach is a powerful alternative to the conventionally used Monte Carlo simulations based on the Boltzmann sampling, and is ideally suited for investigating β-ensembles. Copyright © EPLA, 2013.