The existence and uniqueness of Gevrey regularity solution for a class of nonlinear bistable gradient flows, with the energy decomposed into purely convex and concave parts, such as epitaxial thin film growth and square phase field crystal models, are discussed in this talk. The polynomial pattern of the nonlinear terms in the chemical potential enables one to derive a local in time solution with Gevrey regularity, with the existence time interval length dependent on certain functional norms of the initial data. Moreover, a detailed Sobolev estimate for the gradient equations results in a uniform in time bound, which in turn establishes a global in time solution with Gevrey regularity.