Let $X^N$ be a family of $N\times N$ independent GUE random matrices, $Z^N$ a family of deterministic matrices, $P$ a self-adjoint non-commutative polynomial, that is for any $N$, $P(X^N)$ is self-adjoint, $f$ a smooth function. We prove that for any $k$, if $f$ is smooth enough, there exist deterministic constants $α_i^P(f,Z^N)$ such that $$ \mathbb{E}\left[\frac{1}{N}\text{Tr}\left( f(P(X^N,Z^N)) \right)\right]\ =\ \sum_{i=0}^k \frac{α_i^P(f,Z^N)}{N^{2i}}\ +\ \mathcal{O}(N^{-2k-2}) .$$ Besides the constants $α_i^P(f,Z^N)$ are built explicitly with the help of free probability. In particular, if $x$ is a free semicircular system, then when the support of $f$ and the spectrum of $P(x,Z^N)$ are disjoint, for all $i$, $α_i^P(f,Z^N)=0$. As a corollary, we prove that given $α<1/2$, for $N$ large enough, every eigenvalue of $P(X^N,Z^N)$ is $N^{-α}$-close from the spectrum of $P(x,Z^N)$.