Continental intraplate seismicity seems often episodic, clustered, and migrating. The observed seismicity shows both spatial clustering in seismic zones and scattering across large plate interiors, temporal clustering followed by long periods of quiescence, and migration of seismicity from one seismic zone to another. Here, we explore the complex spatiotemporal patterns of intraplate seismicity using a 3D visco-elasto-plastic finite-element model. The model simulates tectonic loading, crustal failure in earthquakes, and coseismic and postseismic stress evolution. For a laterally homogeneous lithosphere with randomly prespecified perturbations of crustal strength, the model predicts various spatiotemporal patterns of seismicity at different timescales: spatial clustering in narrow belts and scattering across large regions over hundreds of years, connected seismic belts over thousands of years, and widely scattered seismicity over tens of thousands of years. The orientation of seismic belts coincides with the optimal failure directions associated with the assumed tectonic loading. Stress triggering and migration cause spatiotemporal clustering of earthquakes. When weak zones are included in the model the predicted seismicity initiates within the weak zones but then extends far beyond them. If a fault zone is weakened following a large earthquake, repeated large earthquakes can occur on the same fault zone even in the absence of strong tectonic loading. These complex spatiotemporal patterns of intraplate seismicity predicted in this simple model suggest that assessment of earthquake hazard based on the limited historic record may be biased toward overestimating the hazard in regions of recent large earthquakes and underestimating the hazard where seismicity has been low during the historic record.