In this study, an extension of the variational model proposed by Vallabhan and Mustafa (1996) for the settlement analysis of an axially loaded pile embedded in a multi-layered soil profile is presented. In this model, the soil profile and the embedded pile are divided into a number of sub-layers according to the actual number of soil layers observed in the field. The displacement shape function of each soil layer is given as the product of an exponential equation along the pile depth and Bessel’s solution in the radial direction. The displacement relationship in each layer can be derived through transformation matrices. One of the major features of this method is that the total number of pile elements is the same as the total number of soil sub-layers. All the field components, such as displacements, stresses, and strains in the soil, can be calculated by closed-form solutions except the only unknown variable, β , which is a variable as expressed by the modified Bessel’s functions of the second kind, and which can be determined by iteration techniques. Comparisons were made with the results of finite element analysis and field pile-loaded tests. Reasonable results were obtained for all cases.