In this study, the capabilities of four parameter-free models for predicting the heterogeneous azeotropic behavior of ternary and quaternary mixtures were evaluated. The present work was basically extended from that of Lee et al. (1996, 1999) on homogeneous azeotrope determination. The mathematical techniques of Newton-Raphson, pseudo-archlength continuation, and the concept of homotopy were used to accomplish the mathematical computations. Four parameter-free models, Schatchard-Hildebrand (1931), UNIFAC (1975), Vetere-NRTL (1994), and Ash-Wilson (1992), were employed for liquid phase activity coefficient estimations. The ternary and quaternary azeotropic mixtures from Horsley's data book (1973) were adopted in this study. The UNIFAC model has the best performance among these four models. The average absolute relative deviations (AARD%) of azeotropic temperatures and compositions estimated by the UNIFAC model were 1.7%, 1.6%, 15.6%, and 23.9% by comparing with Horsley's data for ternary and quaternary mixtures, respectively. It is very interesting to observe that many homogeneous azeotropes in Horsley's data book appeared to be heterogeneous in our calculations.