In this paper, we present a Hessian recovery based linear finite element method to simulate the molecular beam epitaxy growth model with slope selection. For the time discretization, we apply a first-order convex splitting method and secondorder Crank-Nicolson scheme. For the space discretization, we utilize the Hessian recovery operator to approximate second-order derivatives of a C 0 linear finite element function and hence the weak formulation of the fourth-order differential operator can be discretized in the linear finite element space. The energy-decay property of our proposed fully discrete schemes is rigorously proved. The robustness and the optimal-order convergence of the proposed algorithm are numerically verified. In a large spatial domain for a long period, we simulate coarsening dynamics, where 1/3power-law is observed.