Abstract: This paper develops the new boundary element method (BEM) for accurately calculating additional stresses and settlements in non-horizontally layered foundations. The proposed BEM utilizes the fundamental solutions of multilayered solids with transverse isotropy (TI) under point-concentrated loads. Eight-noded isoparametric elements are used to discretize a core region around the traction area whereas 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 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. Numerical verification shows that the present BEM has very high accuracy and computational efficiency. Numerical examples illustrate the effects of anisotropy and non-horizontal layered structures of the foundations on the elastic fields in detail.