Let \(g:\mathbb R\to \mathbb C\) be a \(C^{\infty}\)-function with all derivatives bounded, and let \(\operatorname{Tr}\) denote the trace on \(n\times n\) matrices. The authors prove for a random matrix \(X_n\) from the Gaussian unitary ensemble (GUE) an asymptotic expansion \[ \frac 1n \operatorname{E} \left[\operatorname{Tr} g\left(\frac{X_n}{\sqrt n}\right)\right]= \frac 1 {2\pi} \int_{-2}^2 g(x) \sqrt{4-x^2}dx + \sum_{j=1}^k \frac{\alpha_j(g)}{n^{2j}}+ O(n^{-2k-2}), \] where \(k\) is an arbitrary positive integer. The coefficients \(\alpha_j(g)\) of the expansion are characterized in terms of Chebyshev polynomials. This gives a proof of a result by N. M. Ercolani and K. D. T.-R. McLaughlin [Int. Math. Res. Not. 2003, No. 14, 755–820 (2003; Zbl 1140.82307)]. A similar asymptotic expansion for the covariance \[ \operatorname{Cov}\left( \operatorname{Tr} f\left(\frac{X_n}{\sqrt n} \right), \operatorname{Tr} g\left(\frac{X_n}{\sqrt n} \right)\right) \] is obtained. Special attention is paid to the case where \(g(x)=\frac{1}{\lambda-x}\) and \(f(x)=\frac 1 {\mu-x}\). In this case, the mean and the covariance considered by the authors are related to the one- and the two-dimensional Cauchy transform of the mean distribution of eigenvalues of \(X_n\).