In this paper a practical numerical method for calculating the settlement response of piles embedded in multilayered soils is presented. Finite beam elements are used for the discretization of piles. The soil mass, which is considered as a semi-infinite or finite medium with each layer being elastic, homogeneous and isotropic, is modeled with the boundary element method and the Steinbrenner's approximation. In addition, the critical case in which the elastic modulus is decreasing with depth is appropriately treated by considering the soil as a half space inhomogeneous medium. The present formulation, besides yielding fully adequate results for design purposes when compared with other approaches, presents a more practical problem modeling and data preparation. Results of numerical examples considering single piles and pile groups are compared with those of other formulations in order to demonstrate the efficiency and accuracy of the proposed formulations.