A computational framework is presented for materials science models that come from energy gradient flows that lead to the evolution of structure involving two or more phases. The models are considered in periodic cells and standard Fourier spectral discretization in space is used. Implicit time stepping is used, and the resulting implicit systems are solved iteratively with the preconditioned conjugate gradient method. The dependence of the condition number of the preconditioned system on the size of the time step and the order parameter in the model (that represents the scaled width of transition layers between phases) is investigated. The framework is easily extended to higher order derivative models, higher dimensional settings, and vector problems. Several examples of its application are demonstrated, including a sixth order problem in three dimensions. A comparison to time-stepping with operator splitting (into convex and concave parts) is done.