This paper presents an efficient and effective scheme for reforming the original Thomson–Haskell transfer matrix method for the stable modeling of surface waves in a layered half-space. The proposed scheme splits the layer propagator matrix and then determines the interface stiffness instead of directly building the propagator matrix as the original Thomson–Haskell method does. For the Love wave, an explicit solution for the interface stiffness is obtained. The reformulation of the Thomson–Haskell method keeps the simplicity of the original method but naturally excludes its exponential growth terms. Hence, the modified Thomson–Haskell method is computationally stable for high frequency and increasing layer thickness. Moreover, both analytical and numerical results indicate the modified Thomson–Haskell method is the most efficient of all the complex arithmetic algorithms available, including the original Thomson–Haskell method and generated reflection and transmission matrix method. The computational speed of the modified Thomson–Haskell method rises by more than 15% for the Love and Rayleigh waves. Furthermore, a discussion of the effect of frequency and layer thickness on the phase velocity of Gutenberg Earth model indicates the modified Thomson–Haskell method is a powerful candidate for efficiently and effectively simulating surface waves in a layered half-space.