We present a 3D hybrid method which combines the Finite Element Method (FEM) and the Spectral Boundary Integral method (SBIM) to model nonlinear problems in unbounded domains. The flexibility of FEM is used to model the complex, heterogeneous, and nonlinear part -- such as the dynamic rupture along a fault with near fault plasticity -- and the high accuracy and computational efficiency of SBIM is used to simulate the exterior half spaces perfectly truncating all incident waves. The exact truncation allows us to greatly reduce the domain of spatial discretization compared to a traditional FEM approach, leading to considerable savings in computational cost and memory requirements. The coupling of FEM and SBIM is achieved by the exchange of traction and displacement boundary conditions at the computationally defined boundary. The method is suited to implementation on massively parallel computers. We validate the developed method by means of a benchmark problem. Three more complex examples with a low velocity fault zone, low velocity off-fault inclusion, and interaction of multiple faults, respectively, demonstrate the capability of the hybrid scheme in solving problems of very large sizes. Finally, we discuss potential applications of the hybrid method for problems in geophysics and engineering.