On the example of the three-dimensional Ising model, we show that nonperturbative renormalization group equations allow one to obtain very accurate critical exponents. Implementing the order ∂ 4 of the derivative expansion leads to ν = 0.632 and to an anomalous dimension η = 0.033 which is significantly improved compared with lower orders calculations.Many problems in high-energy as well as in statistical physics call for nonperturbative methods. On the one hand, several physical systems are described by field theories in their strong coupling regime so that the usual perturbative techniques become troublesome. They fail either because only the first orders of perturbation are computed and do not suffice, or because, even when high orders are known, standard resummation techniques do not provide converged results. On the other hand, some phenomena such as confinement in QCD or phase transitions induced by topological defects are genuinely nonperturbative.Apart from some methods restricted to specific dimensions or situations, very few nonperturbative techniques are available. During the last years, the Wilson approach 1 to the renormalization group (RG) has been turned into an efficient tool 2,3,4 . This nonperturbative RG can be implemented in very general situations and, in particular, in any dimension, so that it has allowed one to study several issues difficult to tackle within a perturbative framework among which the three-dimensional Gross-Neveu model 5 , frustrated magnets 6 , the randomly dilute Ising model 7 , and the Abelian Higgs model 8 .This method relies on a nonperturbative renormalization of the effective action Γ, i.e. the Gibbs free energy. It consists in building an effective action Γ k at the running scale k by integrating out only fluctuations greater than k. At the scale k = Λ, Λ −1 denoting the spacing of the underlying lattice, Γ k coincides with the Hamiltonian H since no fluctuation has yet been taken into account while, at k = 0, it coincides with the standard effective action Γ since all fluctuations have been integrated out. Thus, Γ k continuously interpolates between the microscopic Hamiltonian H and the free energy Γ. The running effective action Γ k follows an exact equation which controls its evolution with the running scale k 2 :where t = ln(k/Λ) and Γ(2) k [φ] is the second functional derivative of Γ k with respect to the field φ(q). In Eq. (1), R k (q) is an infrared cutoff function which suppresses the propagation of the low-energy modes without affecting the high-energy ones.Although exact, Eq. (1) is a functional partial integrodifferential equation which cannot be solved exactly. To handle it, one has to truncate Γ k . A natural and widely used truncation is the derivative expansion, which consists in expanding Γ k in powers of ∂φ, keeping only the lowest order terms. Physically, this truncation rests on the assumption that the long-distance physics of a given model is well described by the lowest derivative terms, the higher ones corresponding to less relevant operators. U...