A computational system is built for conducting deterministic simulations of threedimensional multiple scattering of elastic waves by spherical inclusions.Based on expansion expression of elastic wave fields in terms of scalar and vector spherical harmonics, analytically exact solutions of single scattering and multiple scattering are obtained, implemented and verified. The verification is done by using continuities of displacement and surface traction at the interface between an inclusion and host medium, energy conservation and published results.The scatterer polymerization methodology is extended to three-dimensional multiple scattering solution. By using this methodology, an assemblage of actual scatterers can be treated as an abstract scatterer. This methodology is verified by using different approaches, with or without scatterer polymerization, to solve a physically the same multiple scattering problem.As an application example, band gap formation process for elastic wave propagation in cubic lattice arrangements of spherical scatterers is observed through a series of numerical simulations. Along the direction of the incident wave, scatterer arrangements are viewed as comprising layers of scatterers, within which scatterers form a square grid. Starting from one layer and by increasing the number of layers, near-field forward wave propagation spectra are computed as the number of layers increases.These simulations also demonstrates that the computational system has the capability to simulate multiple scattering solutions of elastic waves in three-dimension. The scatterer polymerization methodology is extended to three-dimensional multiple scattering solution. By using this methodology, an assemblage of actual scatterers can be treated as an abstract scatterer. This methodology is verified by using different approaches, with or without scatterer polymerization, to solve a physically the same multiple scattering problem.As an application example, band gap formation process for elastic wave propagation in cubic lattice arrangements of spherical scatterers is observed through a series of numerical simulations. Along the direction of the incident wave, scatterer arrangements are viewed as comprising layers of scatterers, within which scatterers form a square grid. Starting from one layer and by increasing the number of layers, near-field forward wave propagation spectra are computed as the number of layers increases.These simulations also demonstrates that the computational system has the capability to simulate multiple scattering solutions of elastic waves in three-dimension.