This study extends the analytical solution of folding from a single-layer system to a multi-layer system composed of layers with varying thicknesses, viscosities and matrices, based on the stress distribution of folding with prescribed amplitude and wavelength. The parameters and factors affecting fold behavior were normalized to dimensionless quantities. In addition to single-layer systems, the influences of individual thickness, spacing and viscosity of multiple layers should be of concern. The range of matrix thickness where group folding occurs was initially identified and found to be smaller than one wavelength. Two-dimensional finite element folding simulations were also conducted to confirm the validity of the developed solutions of multi-layer folding. The 2D numerical simulation reveals when the matrix is thinner than one wavelength, the stress distribution within the matrix becomes linear and group folding occurs, and, if the matrix is thicker than one wavelength, the folding of the multi-layer system in fact is identical to independent folding of each layer alone. By means of stress superposition, this study developed theoretical solutions for multi-layer systems. The solutions developed show that the folding of viscous multi-layer and viscous matrix systems is time independent. Based on the developed solution, the effect of thickness, spacing and viscosity are then presented with the other parameters fixed. Finally, the solutions for elastic multi-layer and elastic matrix systems are also proposed.