We present N-body simulations performed in the framework of MOND. This work is based on a numerical resolution of the modified Poisson equation derived from a Lagrangian theory of MOND. This equation is a nonlinear partial differential equation, so we first developed a code that solves this kind of equations using multigrid techniques. We compared the dynamical behaviour of a typical isolated galaxy in Newton CDM model and MOND model. In this approach, pure stellar disc are considered. For the same value of the Toomre parameter (Q_T), galactic discs in MOND develop a bar instability sooner than in the DM model. In a second phase the MOND bars weaken while the DM bars continue to grow by exchanging angular momentum with the halo. The bar pattern speed evolves quite differently in the two models, this affects the position of resonance like the corotation and the peanut. Then we studied the evolution of several galactic discs representing the Hubble sequence, in both models. These simulations lead to a statistical bar frequency which is closer to observations for the MOND than the DM model.