In the terahertz frequency range, the effective (complex-valued) surface conductivity of atomically thick 2D materials such as graphene has a positive imaginary part that is considerably larger than the real part. This feature allows for the propagation of slowly decaying electromagnetic waves, called surface plasmon-polaritons (SPPs), that are confined near the material interface with wavelengths much shorter than the wavelength of the free-space radiation. SPPs are promising ingredients in the design of novel optical applications promising "subwavelength optics" beyond the diffraction limit. There is a compelling need for controllable numerical schemes which, placed on firm mathematical grounds, can reliably describe SPPs in a variety of geometries. In this talk we present a finite element approach for the simulation of SPP structures on a conducting sheet, excited by a plane-wave or electric Hertzian dipole sources. Corresponding analytical results are briefly discussed as well.