This article investigates effects of axial load distribution on buckling loads and their modes of functionally graded (FG) beams including a shear effect for the first time, since all previous studies focused on constant axial load. The parabolic higher order shear theory of Reddy is proposed to consider the shear effect of both thick and thin beams. The gradation of the material is distributed continuously through thickness direction with generalized power-law function. The in-plane axial compressive load is simulated by continuous functions through the axial direction, which are constant, linear, and parabolic functions. The equilibrium equations are developed by Hamilton’s principles, which yields variable-coefficients differential equations. Modified differential quadrature method (DQM) is employed to discretize and mesh the spatial domain and covert the differential equations and boundary conditions to algebraic equations. Algebraic set of equations is formed as a generalized matrix eigenvalue problem, then solved to obtain buckling loads (eigenvalues) and mode-shapes (eigenvectors). Numerical results are presented to illustrate effects of gradation parameter, load types, slenderness ratio and boundary conditions on buckling loads and mode-shapes of functionally graded beams structures. This model can be applied in designing of space-shuttle, spacecraft, nuclear structure, and naval structure subjected to distributed in-plane axial load.