  We describe a numerical method which allows us to go beyond the classical
approximation for the real-time dynamics of many-body systems by approximating
the many-body Wigner function by the most general Gaussian function with
time-dependent mean and dispersion. On a simple example of a classically
chaotic system with two degrees of freedom we demonstrate that this Gaussian
state approximation is accurate for significantly smaller field strengths and
longer times than the classical one. Applying this approximation to matrix
quantum mechanics, we demonstrate that the quantum Lyapunov exponents are in
general smaller than their classical counterparts, and even seem to vanish
below some temperature. This behavior resembles the finite-temperature phase
transition which was found for this system in Monte-Carlo simulations, and
ensures that the system does not violate the Maldacena-Shenker-Stanford bound
$\lambda_L < 2 \pi T$, while the classical dynamics inevitably breaks the
bound.
