This paper is devoted to the development, software implementation, and quantitative estimation of a numerical algorithm based on the Cartesian grid method for the mathematical modeling of shock wave propagation in domains of complex shape with curvilinear boundaries. A detailed description of an algorithm based on the method of “h-boxes” is given. The efficiency of the algorithm is analyzed on the problems of regular and single Mach reflection of a shock wave from a wedge as well as on the problem of shock wave/cylinder interaction.