Over the past decade or more, there has been a great deal of interest in the mechanical, electronic and optical properties of graphene. The interest in the optical properties of monolayer graphene stems largely from two key features of the electronic band structure: it is gapless and the dispersion is linear near the Dirac points, where the two bands touch. Both of these features can be uniquely accessed using terahertz radiation, because, under the right conditions, it can strongly drive both interband and intraband transitions. In this talk, I will describe the formalism and computational results of our density-matrix approach to modelling the nonlinear terahertz response of monolayer graphene as well as unbiased and biased bilayer graphene. As I will show, although the linear dispersion of monolayer graphene does result in a relatively strong nonlinearity due to intraband transitions, in undoped graphene, it is the interplay between the intraband and interband transitions that dominates the nonlinear response. Consequently, we find that the nonlinear terahertz response of bilayer graphene can be greater than that of monolayer graphene, even though it’s dispersion is parabolic.