An incremental variational equation for non-linear motions of multilayer beams composed of n stiff layers and (n - 1) soft cores is derived from the dynamic virtual work equation by an appropriate integration procedure. The kinematical hypotheses of Euler-Bernoulli and Timoshenko beam theories are used to describe the displacement fields of the stiff layers and cores respectively. An efficient solution procedure of incremental harmonic balance method type, with use of finite elements, is developed. To demonstrate its capability, some problems in free non-linear vibrations of multilayer beams are treated by using the procedure. Results are compared with those available in the literature. The effects of damping are also included in this investigation but are described in Part II [1] of this paper in which a number of undamped and damped forced non-linear vibration problems are studied. Results in the form of tables and plots are also presented and comparisons are made with those available in the literature. © 1985.