Reactive distillation with potential liquid-phase split and subsequent formation of a third phase is a highly complex process system nevertheless quite common in the production of useful solvents and biofuels through esterification. The optimal design of such process systems requires the development and solution of reliable and accurate process models that lead to a computationally demanding mathematical problem. In this work, a nonequilibrium (NEQ) model coupled with the orthogonal collocation on finite elements (OCFE) technique is developed for the simulation and optimal design of three-phase reactive distillation systems. The resulting NEQ/OCFE model combines the predictive accuracy of the NEQ model as well as the model reduction and approximation capabilities of the OCFE formulation. Therefore, an accurate but compact in size and thus easier to solve model that accounts for all the physical phenomena, the interactions among the multiple phases, and the occurring chemical reactions becomes available. The NEQ/OCFE model is enriched with an accurate prediction procedure for the identification and tracking of the phase transition boundaries (point of transition for a single liquid phase to two liquid-phase regimes and vice versa) inside the column for varying operating conditions. The model is validated using experimental results and utilized in the optimal design and the dynamic simulation of a staged reactive distillation column for the production of butyl acetate via the esterification reaction of butanol with acetic acid. The optimal column configuration defined by the number of stages in each column section, the feed strategy (single feed or multiple feed points), the location of the feed stages, and the operating conditions are calculated through a rigorous design optimization procedure for tight production purity specifications. At optimal conditions the specific column appears to have both two- and three-phase regions that are separated by the feed stage.