  We present a methodology to construct efficient high-order in time accurate
numerical schemes for a class of gradient flows with appropriate Lipschitz
continuous nonlinearity. There are several ingredients to the strategy: the
exponential time differencing (ETD), the multi-step (MS) methods, the idea of
stabilization, and the technique of interpolation. They are synthesized to
develop a generic $k^{th}$ order in time efficient linear numerical scheme with
the help of an artificial regularization term of the form
$A\tau^k\frac{\partial}{\partial t}\mathcal{L}^{p(k)}u$ where $\mathcal{L}$ is
the positive definite linear part of the flow, $\tau$ is the uniform time
step-size. The exponent $p(k)$ is determined explicitly by the strength of the
Lipschitz nonlinear term in relation to $\mathcal{L}$ together with the desired
temporal order of accuracy $k$. To validate our theoretical analysis, the thin
film epitaxial growth without slope selection model is examined with a
fourth-order ETD-MS discretization in time and Fourier pseudo-spectral in space
discretization. Our numerical results on convergence and energy stability are
in accordance with our theoretical results.
