Limit cycles of perturbed global isochronous center

We apply the averaging method of first order to study the maximum number of limit cycles of the ordinary differential systems of the form ¨x + x = ε (f1(x, y)y + f2 (x, y)) , ¨y + y = ε (g1(x, y)x + g2 (x, y)) , where f1(x, y) and g1(x, y) are real cubic polynomials; f2(x, y) and g2(x, y) are real quadratic polynomials. Furthermore ε is a small parameter.