A new procedure introduced as variational iterative stiffness method (VISM) to overcome the difficulties of the application of variational iteration method to free vibration analysis of higher order beams consist of sixth order differential equation with complex boundaries. A higher order trigonometric shear deformation theory denoted Touratier beam which is considering rotary inertia is used for the mathematical model to perform free vibration analysis of beams with a uniform cross-section, subjected to a compressive axial load, resting on Winkler and Pasternak foundations with various boundary conditions, in this study. Circular frequencies obtained by using VISM are compared with DSM and analytical method (ANM) results based on Touratier and Timoshenko beam theories where a very good agreement is observed. In addition to the obtained circular frequency values, normalised lateral displacements, bending rotation, rotation of the normal to the axis, shear forces and bending moments diagrams of first modes are represented in graphics.