Recently, the molecular beam epitaxy (MBE) has become an important technique for the growth of thin films. By using the MBE technique, it is possible to grow high-quality crystalline materials and form structures with high precision in the vertical direction. There has been a large amount of research interest in the dynamics of the MBE growth. A sufficiently long integration time is necessary in order to detect the epitaxy growth behaviors and to reach the physical scaling regime. To carry out numerical simulations with large time and large computational domain, highly stable and accurate numerical schemes are required. In this talk, our recent works on numerical simulations of these MBE models will be presented.