The main goal of this work is to develop a novel boundary element method (BEM) model for analyzing ultrasonic wave propaga-tion in three-temperature anisotropic viscoelastic porous media. Because of strong nonlinearity of ultrasonic wave propagation in three-temperature porous media problems, the analytical or numerical solutions to the problems under consideration are always challenging, necessitating the development of new computational techniques. As a result, we use a new BEM model to solve such problems. A time stepping procedure based on the linear multistep method is obtained after solving the discretized boundary in-tegral equation with the quadrature rule. The calculation of a double integral is required to obtain fundamental solutions, but this increases the total BEM computation time. Our proposed BEM technique is used to solve the current problem and improve the formulation efficiency. The numerical results are graphed to demonstrate the effects of viscosity and anisotropy on the nonlinear ultrasonic stress waves in three-temperature porous media. The validity, accuracy, and efficiency in the proposed methodology were demonstrated by comparing obtained results to the corresponding solution obtained from finite difference method (FDM).