A numerical method for reaction–diffusion equations with analytic nonlinearity is presented, for which a rigorous backward error analysis is possible. We construct a modified equation, which describes the behavior of the full discretization scheme up to exponentially small errors in the step size. In the construction the numerical scheme is first exactly embedded into a nonautonomous equation. This equation is then averaged with only exponentially small remainder terms. The long-time behavior near hyperbolic equilibria, the persistence of homoclinic orbits and regularity properties are analyzed.