The design of numerical approaches for the molecular beam epitaxy models has always been a hot issue in numerical analysis, in which one of the main challenges for algorithm design is how to establish a high-order time-accurate numerical method with unconditional energy stability. The numerical method developed in this paper is based on the “stabilized-Invariant Energy Quadratization” (S-IEQ) approach. Its novelty is that by adding a very simple linear stabilization term, the difficulty that the original energy potential for the no-slope selection case is not bounded from below can be easily overcome. Then by using the standard format of the IEQ method, we can easily obtain a linear, unconditionally energy stable, and second-order time accurate scheme for solving the system. We further implement various numerical examples to demonstrate the stability and accuracy of the proposed scheme.