The aim of this paper is to provide a mathematical and numerical analysis for a FEM-BEM coupling approximation of a magnetically nonlinear eddy current formulation by using FEM only on the conducting domain, and by imposing the integral conditions on its boundary. The nonlinear relationship between flux density and the magnetic field intensity is given by a physical parameter called magnetic reluctivity, which is assumed to depend on the Euclidean norm of the magnetic induction in the conducting domain. We use the nonlinear monotone operator theory for parabolic equations to show that the continuous formulation obtained for the coupling is a well-posed problem. Furthermore, we use Nédélec edge elements, standard nodal finite elements, and a backward-Euler time scheme, to obtain a fully discrete formulation and to prove quasi-optimal error estimates.