  Over decades, the time evolution of Wigner functions along classical
Hamiltonian flows has been used for approximating key signatures of molecular
quantum systems. Such approximations are for example the Wigner phase space
method, the linearized semiclassical initial value representation, or the
statistical quasiclassical method. The mathematical backbone of these
approximations is Egorov's theorem. In this paper, we reformulate the
well-known second order correction to Egorov's theorem as a system of ordinary
differential equations and derive an algorithm with improved asymptotic
accuracy for the computation of expectation values. For models with easily
evaluated higher order derivatives of the classical Hamiltonian, the new
algorithm's corrections are computationally less expensive than the leading
order Wigner method. Numerical test calculations for a two-dimensional
torsional system confirm the theoretical accuracy and efficiency of the new
method.
