This paper presents a numerical scheme for dynamic analysis of mechanical systems subjected to damping forces that are proportional to fractional derivatives of displacements. These equations appear in the modeling of frequency dependent viscoelastic damping of materials. In the scheme presented, the fractional differential equation governing the dynamics of a system is transformed into a set of differential equations with no fractional derivative terms. Using Laguerre integral formula, this set is converted to a set of first order ordinary differential equations, which are integrated using a numerical scheme to obtain the response of the system. In contrast to other numerical techniques, this method does not require one to store the past history of the response. Numerical studies show that the solution converges as the number of Laguerre node points increase. Further, results obtained using this scheme agree well with those obtained using analytical techniques.