The aim of this paper is to develop a finite element method which allows computingthe buckling coefficients and modes of a non-homogeneous Timoshenko beam.Studying the spectral properties of a non-compact operator,we show that the relevant buckling coefficients correspond to isolatedeigenvalues of finite multiplicity.Optimal order error estimates are proved for the eigenfunctionsas well as a double order of convergence forthe eigenvalues using classical abstract spectral approximation theory for non-compact operators.These estimates are valid independently of the thickness of the beam, whichleads to the conclusion that the method is locking-free.Numerical tests are reported in order to assess the performance of the method.