We study numerically the onset of chaos and thermalization in the Banks-Fischler-Shenker-Susskind (BFSS) matrix model with and without fermions, considering Lyapunov exponents, entanglement generation, and quasinormal ringing. We approximate the real-time dynamics in terms of the most general Gaussian density matrices with parameters which obey self-consistent equations of motion, thus extending the applicability of real-time simulations beyond the classical limit. Initial values of these Gaussian density matrices are optimized to be as close as possible to the thermal equilibrium state of the system. Thus attempting to bridge between the low-energy regime with a calculable holographic description and the classical regime at high energies, we find that quantum corrections to classical dynamics tend to decrease the Lyapunov exponents, which is essential for consistency with the Maldacena-Shenker-Stanford bound at low temperatures. The entanglement entropy is found to exhibit an expected "scrambling" behavior-rapid initial growth followed by saturation. At least at high temperatures the entanglement saturation time appears to be governed by classical Lyapunov exponents. Decay of quasinormal modes is found to be characterized by the shortest timescale of all. We also find that while the bosonic matrix model becomes nonchaotic in the low-temperature regime, for the full BFSS model with fermions the leading Lyapunov exponent, entanglement saturation time, and decay rate of quasinormal modes all remain finite and nonzero down to the lowest temperatures.