The Asymptotic Numerical Method (ANM), (Potier-Ferry et al. [1]) is a robust continuation method for solving non linear fluid or solid mechanical problems that depends on a parameter. Recently, Fan Xu [2], applied the ANM to simulate the non linear behavior film/substrate systems under various kind of applied forces using thin shell elements to take into account the stiff film. Following his pioneer work, in the case of 3D film/system under lateral compressive forces (at the film level), we decided to use the Finite Element Method (FEM) both for the stiff film and the compliant substrate, which allows accounting for various geometries, material properties and boundary conditions, and is more efficient in the case of short instability wavelength. In this paper, we will show that the FreeFem++ Finite Element platform [3] can be successfully used to implement the ANM to simulate wrinkling behavior of 3D film/substrate systems. Wrinkling phenomena in film/substrate systems often requires many degrees of freedom especially for 3D numerical models. Much effort has been spent to develop an efficient MPI parallel code in order to use High Performance Computing capabilities. Numerical results (bifurcation curves, deformed shape, …) will be shown.