Abstract: A new boundary element method (BEM) is developed for accurately calculating the additional stresses and settlements in non-horizontally layered foundations. The proposed BEM utilizes the fundamental solutions for multilayered solids with transverse isotropy (TI) under point-concentrated loads. The eight-noded isoparametric elements are used to discretize a core region around the traction area, whereas the six-noded infinite elements are used to discretize a far-field region beyond the traction area. To avoid calculating the coefficient of the free term for the source point at the strike line between the boundary and the material interface, the discontinuous isoparametric and infinite elements with nodes located within elements are positioned near the strike lines. Numerical methods for non-singular and various singular integrals in the discretized boundary integral equations are developed. The numerical verification shows that the proposed BEM has very high accuracy and computational efficiency. The numerical examples illustrate the effects of anisotropy and non-horizontal layered structures of the foundations on the elastic fields.