Fretting is a major cause of surface damage and fretting fatigue crack initiation at the interface between contact materials subjected to small oscillatory movements. In the present paper, a multi-layered model is developed to analyze fretting fatigue of functionally graded materials (FGMs) with arbitrarily varying shear modulus under plane strain-state deformation. Based on the fact that an arbitrary curve can be approached by a series of continuous but piecewise linear curves, the FGM is divided into several sub-layers and in each sub-layers the shear modulus is assumed to be a linear function while the Poisson’s ratio is assumed to be a constant. With the model, the problem of fretting contact of two similar functionally graded coated cylinders is investigated. By using the transfer matrix method and Fourier integral transform technique, the problem is reduced to two uncoupled Cauchy singular integral equations. The tangential contact pressures in the slip and stick zones are calculated by solving the equations numerically. The results show that appropriate gradual variation of the shear modulus can significantly alter the pressures in the contact zone. This may lead to suppression of fretting fatigue cracks at the edges of the contact zone and thus modify the fretting contact damage.