We present a high-order finite difference method for earthquake cycle simulations within complex geometries that incorporates full dynamic effects and is provably stable. The method is developed for the two-dimensional anti- plane strain problem where a rate-and-state frictional fault is embedded in a heterogeneous elastic half-space. To overcome challenges imposed by the large range of spatial and temporal scales, we use the fault slip rate as a threshold for switching between two different numerical solvers. During the interseismic phase, we use a quasi-static approximation where large linear systems of equations must be solved to allow for large time steps. Previous approaches to simulating the co-seismic phase enforced a rate-and-state friction law that resulted in a stiff system of equations that required a special, semi-implicit time-stepping method. Here, we apply a newly developed non-stiff method that is compatible with the interseismic method and allows for traditional explicit time stepping methods. Our simulations are verified by rigorous convergence tests and code- comparison exercises, and are an important step towards advancing earthquake cycle simulations incorporating more realistic physics.