In this research, an improved first order formulation is presented to study the critical buckling load in functionally graded beams. The formulation has five independent variables in comparison with the Timoshenko theory that has three. The Trefftz criterion is utilized with incremental and fundamental states to define the stability analysis. Virtual work statements are derived for the finite element model where the field variables are interpolated by Lagrange polynomials. The numerical results are compared and verified with other formulations found in literature. Parametric studies are also carried out for buckling behavior due to different slenderness ratios, power-law indices and boundary conditions. Applications of the model to functionally graded materials show the validity of the present approach.