We consider expectations of the form $E [tr h_1(X_1^N)... tr h_r(X_r^N)]$, where $X_i^N$ are self-adjoint polynomials in various independent classical random matrices and $h_i$ are smooth test function and obtain a large $N$ expansion of these quantities, building on the framework of polynomial approximation and Bernstein-type inequalities recently developed by Chen, Garza-Vargas, Tropp, and van Handel. As applications of the above, we prove the higher-order asymptotic vanishing of cumulants for smooth linear statistics, establish a Central Limit Theorem, and demonstrate the existence of formal asymptotic expansions for the free energy and observables of matrix integrals with smooth potentials.