The Parallel Volume Integral Equation Method (PVIEM) is applied for the analysis of elastic wave scattering problems in an unbounded solid containing various types of multiple multilayered anisotropic inclusions. It should be noted that this numerical method does not require the use of the Green's function for the anisotropic inclusion - only the Green's function for the unbounded isotropic matrix is needed. This method can also be applied to solve general elastodynamic problems involving arbitrary shape and number of inhomogeneous and/or anisotropic inclusions. A detailed analysis of the SH wave scattering problem is presented for various types of multiple multilayered orthotropic inclusions. Numerical results are presented for the displacement and stress fields at the interfaces for square and hexagonal packing arrays of various types of multilayered orthotropic inclusions in a broad frequency range of practical interest. The standard parallel programming, such as MPI (message passing interface), was used to speed up computation in the volume integral equation method (VIEM). Parallel volume integral equation method 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 various types of multiple multilayered anisotropic inclusions on displacements and stresses at the interfaces of the inclusions.