Let $U^N$ be a family of $N\times N$ independent Haar unitary random matrices and their adjoints, $Z^N$ a family of deterministic matrices, and $P$ a self-adjoint noncommutative polynomial, i.e. for any $N$, $P(U^N,Z^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(U^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. As a corollary, we prove that given $α<1/2$, for $N$ large enough, every eigenvalue of $P(U^N,Z^N)$ is $N^{-α}$-close to the spectrum of $P(u,Z^N)$ where $u$ is a $d$-tuple of free Haar unitaries. We also prove the convergence of the norm of any polynomial $P(U^N\otimes I_M, I_N\otimes Y^M)$ as long as the family $Y^M$ converges strongly and that $M\ll N \ln^{-3}(N)$.