A new numerical model for analyzing the buckling of columns with variably distributed axial loads is proposed in this paper. The presented method transforms the traditional eigenvalue problem into an initial boundary value problem which can be solved by numerical integrations. The buckling load is determined by a proposed two-step iterative procedure. Because of its concise and easy-to-implement form, the proposed method can be coded and used for analysis with a minimum requirement for the in-core memory storage of a computer. The proposed model is verified for a column with uniformly distributed axial loads, and the application of the method is demonstrated by analyzing columns with variably distributed axial loads and columns with varying cross-sections. Some interesting conclusions are made.