The plasmonics of graphene and other two-dimensional materials has attracted enormous amounts of attention in the scientific literature over the past decade. Both the possibility of exciting plasmons in the terahertz to midinfrared regime, and the active tunability of graphene via electrical gating or chemical doping has generated a great deal of excitement among engineers seeking sensing devices which operate in this regime. Consequently, there is significant demand for robust and highly accurate computational capabilities which can incorporate such materials. Standard volumetric approaches can answer this demand, but require vast computational resources in exchange. Here we describe an algorithm which addresses this issue in two ways, first, we model the graphene layer with a surface current which is applicable to a wide class of two-dimensional materials. In addition, we reformulate the governing volumetric equations in terms of surface quantities using Dirichlet-Neumann Operators. These surface equations can be numerically simulated in an efficient, stable, and accurate fashion using a novel High-Order Perturbation of Envelopes methodology. We utilize an implementation of this algorithm to study absorbance spectra of TM polarized plane-waves scattered by a periodic grid of graphene ribbons.