Taking the extended Mindlin solution of a vertical point load acting in the interior of a multilayered elastic medium as the basic solution, an indirect boundary element method (BEM) was adopted to analyze pile foundations embedded in a multilayered elastic soil. The numerical analysis was carried out by the compiled program for a pile foundation in an homogeneous soil and a pile foundation in a multilayered soil. Reasonably good agreement is obtained between the present and existing solutions. Compared with the FEM software ABAQUS, the computational speed of BEM is much faster than that of FEM. Applying BEM in practice is convenient. It is pointed out that the exact measurement of elastic modulus is the key factor for applying BEM in geotechnical engineering.