Rolling of thin sheets generally induces flatness defects due to thermo-elastic deformation of rolls. This leads to heterogeneous plastic deformations throughout the strip width and then to out of plane displacements to relax residual stresses. In this work we present a new numerical technique to model the buckling phenomena under residual stresses induced by rolling process. This technique consists in coupling two finite element models: the first one consists in a three dimensional model based on 8-node tri-linear hexahedron which is used to model the three dimensional behaviour of the sheet in the roll bite; we introduce in this model, residual stresses from a full simulation of rolling (a plane-strain elastoplastic finite element model) or from an analytical profile. The second model is based on a shell formulation well adapted to large displacements and rotations; it will be used to compute buckling of the strip out of the roll bite. We propose to couple these two models by using Arlequin method. The originality of the proposed algorithm is that in the context of Arlequin method, the coupling area varies during the rolling process. Furthermore we use the asymptotic numerical method (ANM) to perform the buckling computations taking into account geometrical nonlinearities in the shell model. This technique allows one to solve nonlinear problems using high order algorithms well adapted to problems in the presence of instabilities. The proposed algorithm is applied to some rolling cases where ``edges-waves'' and ``center-waves'' defects of the sheet are observed.