In this article, a numerical solution methodology is presented to study the postbuckling configurations and free vibrations of Timoshenko beams undergoing postbuckling. The effect of geometrical imperfection is taken into account, and the analysis is carried out for different types of boundary conditions. Based on Hamilton's principle, the governing equations and corresponding boundary conditions are derived. After introducing a set of differential matrix operators that is used to discretize the governing equations and boundary conditions, the pseudo-arc length continuation method is applied to solve the postbuckling problem. Then, the problem of free vibration around the buckled configurations is solved as an eigenvalue problem using the solution obtained from the nonlinear problem in the previous step. This study shows that, when the axial load in the postbuckling domain increases, the vibration mode shape of buckled beam corresponding to the fundamental frequency may change. Another finding that can be of great technical interest is that, for all types of boundary conditions and in both prebuckling and postbuckling domains, the natural frequency of imperfect beam is higher than that of ideal beam. Also, it is observed that, by increasing the axial load, the natural frequency of both ideal and imperfect beams decreases in the prebuckling domain, while it increases in the postbuckling domain. The reduction of natural frequency in the transition area from the prebuckling domain to the postbuckling domain is due to the severe instability of the structure under the axial load.