We compute the critical exponents ν, η and ω 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 O(∂4)]. We analyze the behavior of this approximation scheme at successive orders and observe an apparent convergence with a small parameter, typically between 19 and 14, 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 long-standing controversy between Monte Carlo estimates and experiments for the specific heat exponent α, our results are compatible with those of Monte Carlo but clearly exclude experimental values.