This paper nonlinearly models cantilever-based functionally graded magneto-electro-elastic energy harvesters (FGMEEEH) for the first time. The coupled magneto-electro-mechanical model is obtained on the basis of the Euler-Bernoulli beam theory. A hybrid procedure including Ritz's method is then utilized to generate reduced order models for both asymmetric unimorph and symmetric bimorph configurations. The resulting sets of initial value problems, whose convergence will be examined, are then analytically solved using the method of multiple time scales for both the free vibrations and primary resonance cases. The analytical time-histories of the system are compared by those obtained numerically and excellent agreements between them are observed. In addition, simplifying the frequency response function of the system in the primary resonance case, the present findings are validated by those available in the literature for linear unimorph systems. The influences of the harvester configurations, base acceleration amplitude, the value of the tip mass, the material gradation index as well as the resistances of the piezoelectric and magnetic circuits on the nonlinear response of the system are also studied in detail. It is observed that FGMEEEEHs enjoy much more efficiency in comparison to piezoelectric-based systems.