The parallel volume integral equation method (PVIEM) is applied for the analysis of elastic wave scattering problems in an unbounded isotropic solid containing multiple multilayered anisotropic elliptical inclusions. This recently developed numerical method does not require the use of Green’s function for the multilayered anisotropic inclusions; only Green’s function for the unbounded isotropic matrix is needed. This method can also be applied to solve general two- and three-dimensional elastodynamic problems involving inhomogeneous and/or multilayered anisotropic inclusions whose shape and number are arbitrary. A detailed analysis of the SH wave scattering is presented for multiple triple-layered orthotropic elliptical inclusions. Numerical results are presented for the displacement fields at the interfaces for square and hexagonal packing arrays of triple-layered elliptical inclusions in a broad frequency range of practical interest. It is necessary to use standard parallel programming, such as MPI (message passing interface), to speed up computation in the volume integral equation method (VIEM). Parallel volume integral equation method as a pioneer of numerical analysis enables us to investigate the effects of single/multiple scattering, fiber packing type, fiber volume fraction, single/multiple layer(s), multilayer’s shape and geometry, isotropy/anisotropy, and softness/hardness of the multiple multilayered anisotropic elliptical inclusions on displacements at the interfaces of the inclusions.