### Literature Overview
In his model, both the cross-section bending stiffness and the axial load variation are described through polynomial expressions, and the stiffness matrix for the member includes the effects of the axial loading [16]. The differential equation, which governs the phenomena, can be transformed into a system of algebraic equations with unknown coefficients. The buckling load is found as the load that makes the determinant of the stiffness matrix equal to zero [16]. The Variational Iteration Method (VIM) is a powerful method in the analysis of various engineering problems, and for solving nonlinear ordinary and partial differential equations and integral equations [17]. [PERSON] and [PERSON] applied VIM for the determination of critical buckling loads for Euler columns with constant and variable cross-sections. The study assumes various buckling cases with different boundary conditions and with different variations of cross-sections [17]. For the VIM exemplification method, there are presented a constant cross-section with an initial approximation of the linear ordinary differential equation is a cubic polynomial function with unknown coefficients. The variation of the cross-section is due to the variation of the moment of inertia of the column along the length of the Euler-type element. The material is assumed to be homogeneous therefore the variation of flexural rigidity depends only on the variation of the cross-section. There are proposed two variations of the flexural rigidity, thus exponential variation and by a power low, admitting a linear, a quadratic, and a cubic variation of the cross-section. The results are compared with the exact analytical solutions [17]. [PERSON] et al. used the variational method to define the static problem, in which the axial displacement field is to be computed [9]. [PERSON] [18] studied the elastic buckling of steel columns under axial compression with square, rectangle, and circular cross-sections, and different boundary conditions, i.e., fixed-free and pinned-pinned boundary conditions. For the solution along numerical computations, finite element modelling (FEM) has been employed. In this study, an analysis has been made regarding the effects of the boundary conditions, cross-sections, and slenderness ratios on the buckling loads of the steel columns [18]. Advanced numerical analysis can be performed with ANSYS Parametric Design Language codes, in order to determine the critical buckling loads. [PERSON] and [PERSON] investigated circular variable cross-section concrete columns with a triangular variation of the cross-section admitting different ratios for the prismatic and the tronconic segments along the height of the columns. The analysis confirms that significant material economy and critical buckling load maximization can be obtained [12]. [PERSON] et al. [11] proposed a theoretical elastic stability analysis of columns with variable cross-sections subjected to combined distributed and concentrated axial load. The minimization of the total potential energy technique is employed. The variation of the moment of inertia and the load intensity is governed by a power law of the distance along the column length. A simplified concept of the modified slenderness ratio is used to compute the critical stress for columns. The minimum potential energy principle is used to establish the necessary equations, which take into consideration the first two terms of the Fourier series for the deflection shape of the columns, which satisfies the geometric conditions [11]. [PERSON] suggested that one should not disregard the variation of non-dimensional critical buckling load with the variation of temperature field, which is a current research topic [19]. A localized differential quadrature method (LDQM) has a high potential to implement in the buckling analysis of axially graded non-uniform columns with elastic restraints. LDQM can be implemented for solving generalized eigenvalue problems with high accuracy and less computational effort [10].