Understanding the mechanics and rheological properties of the solid Earth requires integrating observations of earthquake cycle deformation with advanced numerical models of the underlying processes. In this article, we introduce a new numerical modeling tool designed for two-dimensional subduction zone geometries that solves the equations for mechanical equilibrium in the lithosphere-asthenosphere system to predict displacements, strains, and stresses throughout the medium. Using a Boundary Element Method (BEM) framework, we reduce the governing partial differential equations to a set of coupled ordinary differential equations, simulating periodic earthquake cycles in a viscoelastic medium with spatially variable power-law viscous rheologies and rate-dependent fault friction. Our approach overcomes key limitations of existing BEM models by ensuring (1) mechanical consistency across timescales, from individual earthquake cycles to geological periods, and (2) precise stress transfer calculations in a viscoelastic medium. We also demonstrate that for spatially heterogeneous linear viscoelastic materials, the coupled ordinary differential equations can be simplified to an eigenvalue problem, significantly enhancing computational efficiency. These advancements offer a powerful tool for predicting spatiotemporal patterns of surface displacement given complex mantle structures and lay the foundation for high-dimensional inverse problems to infer constitutive properties of the lithosphere-asthenosphere system.