© Int. J. of GEOMATE. All rights reserved. An iterative algorithm to construct a displacement based finite element method for analyzing axially loaded pile embedded in finite depth of elasto-plastic soil is presented. The investigation herein is conducted on the condition of shape function by which exact solutions may be reproduced at the nodal points regarding to a few number of elements. The examined shape functions which satisfy the homogeneous governing equations in elastic and plastic soil through bisection iterative algorithm are introduced to obtain the so-called exact element stiffness matrix via total potential energy principle. Numerical examples of elastostatic pile embedded in elasto-plastic Winkler foundation illustrate the accuracy of proposed element compare with conventional finite element shape functions. Axial force and displacement solutions show very good agreement with data from the available literature i.e. the exact nodal displacement solution is obtained correspond to any point load level even with a single element mesh employed.