An accurate and efficient time-dependent wave packet method is proposed for solving the product state-resolved reaction probabilities of the tetratomic reactive system. In this method, the entire scattering process is divided into the interaction region and multiple asymptotic regions, sharing the same spirit as the interaction-asymptotic region decomposition (IARD) approach in a triatomic reactive scattering process. The hyperspherical coordinate is adopted in the interaction region, while the corresponding Jacobi coordinate is employed in each asymptotic region. Therefore, in this IARD method, the "coordinate problem", the difficulty of expressing the wave function in the entire region using a single coordinate system, can be effectively avoided, and only a very small number of the grid points (or the basis functions) are required. For the numerical illustration, the typical tetratomic reaction H2 + OH with zero total angular momentum is calculated, and compared with other quantum wave packet methods. Our proposed IARD method for the tetratomic reactive system is much more efficient and accurate.