Compact and high speed electromechanical systems lead to higher and higher levels of multiaxial mechanical stress, that may strongly change the magnetic behavior of materials, making the development of highly accurate magnetic models a very important task. Among all available magnetoelastic models, multiscale approaches seem to be the most promising.The coupling effect is introduced at the single crystal scale, mainly attributed to the evolution of magnetic domain structures under magneto-mechanical loading. All these models use however a formulation of free energy where the energetic term describing magneto-elastic coupling is linearly stress dependent, that hardly allow (using artificial mechanisms) providing non-monotonic stress effect on the magnetic behavior. The proposition detailed in this paper is to consider a second order stress term in the free energy expression, allowing a linear dependance of magnetostriction tensor with stress to be defined and providing the non-monotonous stress effect. This introduction leads to a more complex description of magnetoelastic effect that needs the identification of a large number of complementary material constants. In this paper, developments are made in the frame of cubic symmetry for first order magneto-elastic term (joining the classical description) and considering an isotropic second order stress effect for the sick of simplicity. This simplification leads to only two additive physical constants to be identified. An identification procedure is proposed and applied to model the magnetoelastic behavior a non-oriented (NO) 3wt%silicon-iron electrical steel.