In this paper, a numerical method for solution of the free vibration of beams governed by a set of second-order ordinary differential equations of variable coefficients, with arbitrary boundary conditions, is presented. The method is based on numerical integration rather than the numerical differentiation since the highest derivatives of governing functions are chosen as the basic unknown quantities. The kernelsof integral equations turn out to be Green's function of corresponding equation with homogeneous boundary conditions. The accuracy of the proposed method is demonstrated by comparing the calculated results with those available in the literature. It is shown that good accuracy can be obtained even with a relatively small number of nodes.