High-order thin layer method for viscoelastic wave propagation in stratified media. Wave propagation in stratified media has wide applications in petroleum exploration, geophysical inversion, nondestructive evaluation of highway and airport pavement structures, countermine technology, structural health monitoring, and vehicle weigh-in-motion systems. A high-order thin layer method is developed to improve the accuracy and robustness of thin layer method for analyzing viscoelastic wave propagation in stratified media. It approximates the stiffness matrix involving transcendental functions by truncating the Taylor series of the stiffness matrix to the fourth order term. A generalized eigenvalue problem is then formulated, which allows an efficient numerical algorithm to be designed in a computer program DynaThinLayer. The high-order thin layer method is compared to the thin layer method, the stiffness matrix method and the exact Zoeppritz method in terms of accuracy and efficiency. It is shown that the high-order method requires the same level of computation time as does the thin layer method, but the former produce more accurate result than the latter. The computation time of the high-order thin layer method increases linearly when the number of layers increases. The new method is most applicable to situations where a large number of layers is involved or to situations where some natural layers have large thickness.