Various computational methods have been studied with respect to their suitability for obtaining very accurate solutions of plate vibration problems, especially for the higher modes. Because of the interest in the higher modes, also higher order effects such as transverse shear deformation and rotational inertia are considered. The Rayleigh-Ritz method with global trial functions appeared to be a suitable choice. To reach a high convergence rate in order to obtain accurate solutions, the complementary boundary conditions formulated by Baruh and Tadikonda should be satisfied. This can be accomplished when polynomials are used as trial functions. When the polynomials are not properly chosen, the algorithm is not numerically stable. It is shown that orthogonalization of the polynomials by means of the Gram-Schmidt process results in a numerical stable process. For free-free boundary conditions, these orthogonal polynomials are the well known Legendre polynomials. For other boundary conditions the resulting polynomials are very similar to the Legendre polynomials. Because of the very high convergence rates, these methods are suitable for obtaining accurate solutions. The numerical stability guarantees that also the higher modes can be calculated.