  Many-body chaos has emerged as a powerful framework for understanding
thermalization in strongly interacting quantum systems. While recent analytic
advances have sharpened our intuition for many-body chaos in certain large $N$
theories, it has proven challenging to develop precise numerical tools capable
of exploring this phenomenon in generic Hamiltonians. To this end, we utilize
massively parallel, matrix-free Krylov subspace methods to calculate dynamical
correlators in the Sachdev-Ye-Kitaev (SYK) model for up to $N = 60$ Majorana
fermions. We begin by showing that numerical results for two-point correlation
functions agree at high temperatures with dynamical mean field solutions, while
at low temperatures finite-size corrections are quantitatively reproduced by
the exactly solvable dynamics of near extremal black holes. Motivated by these
results, we develop a novel finite-size rescaling procedure for analyzing the
growth of out-of-time-order correlators (OTOCs). We verify that this procedure
accurately determines the Lyapunov exponent, $\lambda$, across a wide range in
temperatures, including in the regime where $\lambda$ approaches the universal
bound, $\lambda = 2\pi/\beta$.
