  We compute the critical exponents $\nu$, $\eta$ and $\omega$ of $O(N)$ models
for various values of $N$ by implementing the derivative expansion of the
nonperturbative renormalization group up to next-to-next-to-leading order
[usually denoted $\mathcal{O}(\partial^4)$]. We analyze the behavior of this
approximation scheme at successive orders and observe an apparent convergence
with a small parameter -- typically between $1/9$ and $1/4$ -- compatible with
previous studies in the Ising case. This allows us to give well-grounded error
bars. We obtain a determination of critical exponents with a precision which is
similar or better than those obtained by most field theoretical techniques. We
also reach a better precision than Monte-Carlo simulations in some physically
relevant situations. In the $O(2)$ case, where there is a longstanding
controversy between Monte-Carlo estimates and experiments for the specific heat
exponent $\alpha$, our results are compatible with those of Monte-Carlo but
clearly exclude experimental values.
