A recent experimental study [1] has demonstrated the alignment of phase separated domains across hundreds of bilayer units in multicomponent stacked lipid bilayers. The origin of this alignment is the interlamellar coupling of laterally phase separated domains. Here, we develop a theoretical model that presents the energetics description of this phenomenon based on the minimization of the free energy of the system. Specifically, we use solution theory to estimate the competition between energy and entropy in different stacking configurations. The model furnishes an elemental phase diagram, which maps the domain distributions in terms of the strength of the intra- and inter-layer interactions and estimates the value of inter-layer coupling for complete alignment of domains in the stacks of five and ten bilayers. The area fraction occupied by co-existing phases was calculated for the system of the minimum free energy, which showed a good agreement with experimental observations.