Exploring the earthquake cycle for large, complex tectonic boundaries that deform over thousands of years requires the development of sophisticated and efficient models. In this paper we introduce a semianalytic three-dimensional (3-D) linear viscoelastic Maxwell model that is developed in the Fourier domain to exploit the computational advantages of the convolution theorem. A new aspect of this model is an analytic solution for the surface loading of an elastic plate overlying a viscoelastic half-space. When fully implemented, the model simulates (1) interseismic stress accumulation on the upper locked portion of faults, (2) repeated earthquakes on prescribed fault segments, and (3) the viscoelastic response of the asthenosphere beneath the plate following episodic ruptures. We verify both the analytic solution and computer code through a variety of 2-D and 3-D tests and examples. On the basis of the methodology presented here, it is now possible to explore thousands of years of the earthquake cycle along geometrically complex 3-D fault systems.