Two-dimensional materials constitute an exciting platform for nonlinear optics with large nonlinearities that are tunable by gating. Hence, gate-tunable harmonic generation and intensity-dependent refraction have been observed in, e.g., graphene and transition-metal dichalcogenides, whose electronic structures are accurately modeled by the (massive) Dirac equation. We exploit on the simplicity of this model and demonstrate here that arbitrary nonlinear response functions follow from a simple iterative approach. The power of this approach is illustrated by analytical expressions for harmonic generation and intensity-dependent refraction, both computed up to ninth order in the pump field. Moreover, the results allow for arbitrary band gaps and gating potentials. As illustrative applications, we consider (i) the gate dependence of third- and fifth-harmonic generation in gapped and gapless graphene, (ii) the intensity-dependent refractive index of graphene up to ninth order, and (iii) the intensity dependence of high-harmonic generation.