Twisted bilayer graphene (TBG) is obtained by stacking two sheets of graphene on top of each other with a relative twist. At incommensurate twist angles, TBG is not periodic and thus does not admit a Brillouin zone or periodic branches of spectrum. Instead, the atoms form a structure which is approximately periodic with respect to the so-called moire lattice, whose unit-cell area is inversely proportional to the square of the twist angle. Thus, at small twist angles, it becomes intractable to simulate electron dynamics directly from first principles, making it essential to derive accurate effective models that reduce the complexity of the problem. The first-order continuum partial differential equation (PDE) model proposed by Bistritzer and MacDonald accurately describes the single-particle electronic properties of twisted bilayer graphene at small twist angles. Here, higher-order corrections to the Bistritzer–MacDonald (BM) model were obtained via a systematic multiple-scales expansion. It was proven that the solution of the resulting higher-order PDE model accurately approximates the corresponding tight-binding wave function under a natural choice of parameters and given initial conditions that are spectrally localized to the monolayer Dirac points. Numerical simulations of tight-binding and continuum dynamics demonstrate the validity of the higher-order continuum model. Symmetries of the higher-order models are also discussed.