In this paper, a two-dimensional phase field crystal model of graphene and hexagonal boron nitride (hBN) is extended to include out-of-plane deformations in stacked multilayer systems. As proof of principle, the model is shown analytically to reduce to standard models of flexible sheets in the small deformation limit. Applications to strained sheets, dislocation dipoles, and grain boundaries are used to validate the behavior of a single flexible graphene layer. For multilayer systems, parameters are obtained to match existing theoretical density functional theory calculations for graphene/graphene, hBN/hBN, and graphene/hBN bilayers. More precisely, it is shown that the parameters can be chosen to closely match the stacking energies and layer spacing calculated by Zhou et al. [Phys. Rev. B 92, 155438 (2015)PRBMDO1098-012110.1103/PhysRevB.92.155438]. Further validation of the model is presented in a study of rotated graphene bilayers and stacking boundaries. The flexibility of the model is illustrated by simulations that highlight the impact of complex microstructures in one layer on the other layer in a graphene/graphene bilayer.