An efficient calculation method of potential distribution in the presence of an embedded body in multi-layer earth has been proposed by expanding the method of image with a consideration of multiple reflection between the ground surface and each underground boundary. For this method, when solving boundary integral equation with the potential of embedded body surface as only one unknown, i.e., when obtaining discretization equation, ordinary boundary element program developed for analyzing the finite closed region can be used. As an example, numerical calculation was conducted for the two-layer earth. The analysis expression of potential distribution in the case of the certain embedded body in two-layer earth has never published. Accordingly, the calculated results were compared with those by the integral equation method. As a result, it was concluded that the primary potential obtained from the present method agreed well with that obtained from the integral equation method. However, there was a disregarded difference in the secondary potential. For confirming the effectiveness, it was necessary to compare with another numerical calculation method, such as finite element method. 5 refs., 5 figs.