Many natural and engineering problems follow gradient flow structures in the sense that systems evolve to decrease certain energy. The dynamics of most of these gradient systems are complicated and hence numerical methods are called for. There are at least two desirable features for numerical algorithms for gradient flows with long evolution process: efficient higher order in time, and long-time stability. We present a class of efficient higher-order energy stable methods for a class of gradient flows based on the exponential time differencing (ETD) method combined with multi-step methods and interpolation. As a specific example, we present a third order ETD based scheme for thin film epitaxial growth model together with numerical results establishing the convergence and stability of the scheme, and the ability of the scheme to capture long-time scaling properties of the system.