This paper uses constraint satisfaction and optimization to find least energy solutions of a solid, elastic frictional model for parallel folding. Such a model is representative of multilayer geological systems undergoing buckling deformation and modelling the evolution of folds poses a significant problem. Simplifying the model down to a two-layer formulation and, assuming the geometry of the whole layered material is governed by this, the behaviour of the central interface is modelled using a number of points whose movement is constrained. The small-deflection model, with a linear foundation and no apparent lock-up criteria, closely matches the sequential nature seen in experiments. By adding a hardening non-linear foundation stiffness and linearly increasing the overburden pressure, the destabilization and restabilization of the experimental load-deflection plots is clearly observed.