Based on the meshless natural element method and incremental theory of plasticity, a new algorithm for elastoplastic analysis of a symmetric structure is put forward. The meshless method is a recently developed numerical method, which only requires the nodal information and can effectively overcome the difficulties brought about by the meshes. As widely known, conventional meshfree methods suffer from the numerical difficulties in dealing with the essential boundary conditions. Differing from most of meshfree methods, the natural element method (NEM) is a Galerkin-based meshless method that is built upon the notion of the natural neighbor interpolation. It is noteworthy that, there are two natural neighbor interpolants: Natural neighbor-based Sibson interpolation and Laplace interpolation (non-Sibsonian interpolation). The interpolation relies on the concepts of Delaunay triangulations and Dirichlet tessellations of a set of nodes to build the shape function. Compared with the moving least squares approximation widely used in many meshless methods, the natural neighbor interpolation does not involve complex matrix inversion, needs no artificial parameter but can improve the computational efficiency. Because the constructed shape functions from the natural neighbor interpolation possess the Kronecker delta function property, the essential boundary conditions can be imposed accurately without any special techniques. Meanwhile, the NEM does not require extra efforts to generate the background because it utilizes a set of Delaunay triangles, which are automatically identified in process of the basis function definition, as its background cell. Therefore, the NEM is a most promising numerical method which can be widely applied to nonlinear analysis of materials, three-dimensional calculation and large deformation problems. Generally, plasticity would make many projects more reasonable, and would fully utilize intensity potentially of structures. The NEM has been successfully used to solve elastoplastic problems of plane structures. In the present work, the NEM is applied to elastoplastic analysis of axisymmetric structures. Firstly, a brief description of the rudiments of Voronoi diagrams and Delaunay triangles in the context of natural neighbor interpolation is delineated. The basic principle of Sibson natural neighbor interpolation approach is also presented. A three-dimensional axisymmetric structure is created by the rotation of the cross section around the axis of symmetry and a set of properly scattered nodes within the cross section are employed. In the analysis of elastoplastic problems, the increments of stress and strain for axisymmetric problems are utilized to characterize the elastoplastic constitutive relationship. Owing to the nonlinear nature of plastic deformation, the computation at each load step employs an iterative solution method, either the Newton-Raphson method or the modified Newton-Raphson method. In order to reduce the expensive operations of stiffness matrix evaluation and factorization, the modified Newton-Raphson method with initial elastic stiffness matrix is adopted in the paper. The tangent predictor-radial return algorithm is applied to integrate the constitutive equations. The elastoplastic analyses of two classical axisymmetric structures have been successfully carried out and the fact that the present results show a good agreement with the analytical solutions demonstrate the effectiveness and accuracy of the proposed method for axisymmetric elastoplastic analyses.