Abstract: The natural element method (NEM) has some disadvantages such as: the existing interpolation function can only be applied linearly on boundaries, the low precision of stress and strain results, and the least square method is required to generate the final results. Recognizing the above disadvantages, this paper develops a new interpolation method, which uses Shepard shape function (0 level MLS shape function) as basic function. The high-power basic function is constructed through Taylor expansion, thus the continuity of NEM is maintained. This method has features of interpolation through points and interpolation C∝ on boundaries. These merits improve the accuracy of stress and strain results, which can be obtained directly through interpolation at Gauss integral points. It is shown the new method can improve the efficiency and accuracy of calculation significantly.