An accurate and easily extendable method to deal with lattice dynamics of solids is offered. It is based on first-principles molecular dynamics simulations and provides a consistent way to extract the best possible harmonic - or higher order - potential energy surface at finite temperatures. It is designed to work even for strongly anharmonic systems where the traditional quasiharmonic approximation fails. The accuracy and convergence of the method are controlled in a straightforward way. Excellent agreement of the calculated phonon dispersion relations at finite temperature with experimental results for bcc Li and bcc Zr is demonstrated.