In this paper, a high-order cylindrical beam model, where the shear deformation is taken into account, will be used to analyze the buckling behaviors of the functionally graded cylindrical beams with radially and axially varying material inhomogeneities. The coupled governing equations for buckling of a cylindrical beam under axial compression are derived, which can be translated into a single differential equation by introducing an auxiliary function. The shifted Chebyshev polynomials are used to compute the critical buckling loads for kinds of boundary conditions. By comparing with the three-dimensional solutions for buckling of homogeneous circular beams, the validity of the introduced model is confirmed. Two typical material property distributions defined by the exponential- and power-law are considered. A parametric study is carried out to investigate the effects of material gradient indexes on the critical buckling loads.