The kernel functions of the boundary integral equations for thin plate were determined by the fundamental solutions for an infinite thin plate. By the discretization of the plate interior and boundary as well as the assumption of the distribution states of plate nodes and foundation reaction forces, the BEM equations of the plate can be established. Meanwhile, based on the analytical layer element solutions for layered foundations, the flexibility matrix of the foundation was obtained by a two-dimensioned Guass-Legendre quadrature. Taking into account the compatible conditions of the displacements at the soils-plate interface, the global BEM equations for the interaction problem between the layered foundation and the thin plate were then established. The solutions for the problem were further obtained by solving the global BEM equations. The accuracy of the present method was verified by comparing existing solutions with the numerical results obtained from the corresponding FORTRAN program in this study. It is observed from numerical examples that when a square thin plate is placed on a foundation, the settlement difference between the two lines perpendicular to y or x coordinate decreases as they approach the center of the plate, and the difference decreases with the decrease of the plate-soil stiffness ratio. Furthermore, the settlement discrepancy between the plate center and the midpoint of the long side is unapparent with the increasing length-width ratio, and the similar variation trend can be found between the midpoint of the wide side and angular point.