  We report a theoretical study of the non-linear magnetoelectric response of
GdFeO$_3$ through an analytical approach combined with a Heisenberg model which
is fitted against first-principles calculations. Our theory reproduces the
non-linear change of polarization under applied magnetic field reported
experimentally such that it allows to analyze the origin of the large responses
in the different directions. We show that the non-linear character of the
response in these materials originates from the fact that the antiferromagnetic
order of Gd atoms changes non-linearly with respect to the applied magnetic
field. Our model can be generalized to other materials in which the
antiferromagnetic ordering breaks inversion symmetry.
