A general analytic approach is proposed for nonlinear eigenvalue problems governed by nonlinear differential equations with variable coefficients. This approach is based on the homotopy analysis method for strongly nonlinear problems. As an example, a beam with arbitrary variable cross section acted on by a compressive axial load is used to show its validity and effectiveness. This approach provides us with great freedom to transfer the original nonlinear buckling equation with variable coefficients into an infinite number of linear differential equations with constant coefficients that are much easier to solve. More importantly, it provides us with a convenient way to guarantee the convergence of solution series. As an example, the beam displacement and the critical buckling load can be obtained for arbitrary variable cross sections. The influence of nonuniformity of moment of inertia is investigated in detail and the optimal distributions of moment of inertia are studied. It is found that the critical buckling load of a beam with the optimal distribution of moment of inertia can be approximately 21–22% larger than that of a uniform beam with the same average moment of inertia. Mathematically, this approach is rather general and thus can be used to solve many other linear/nonlinear differential equations with variable coefficients.