Isobaric, isenthalpic (PH) flash is challenging for multiphase non-isothermal flow simulation using an equation of state (EOS). This is because the number of equilibrium phases is unknown in temperature and composition space, and because the system of equations in PH flash becomes nearly degenerate for narrow-boiling fluids. The term "narrow-boiling" is used in the literature to refer to enthalpy that is sensitive to temperature. The primary objective of this research is to develop the multiphase PH-flash algorithm integrated with stability analysis that resolves the two technical challenges mentioned above. The secondary objective is to present a new analysis of narrow-boiling behavior by coupling energy and phase behavior equations through the temperature dependency of K values. The thermodynamic model used is the Peng-Robinson EOS with the van der Waals mixing rules. PH flash in this research is formulated by use of the tangent-plane-distance function, in which phase split computation is integrated with phase-stability analysis. The formulated PH flash is solved by the direct-substitution algorithm with an arbitrary number of sampling compositions (N-S), at which phase stability is measured during the iteration. The number of equilibrium phases is not required to be fixed in the new algorithm. Results in case studies show that the new algorithm can robustly handle phase appearance/disappearance with narrow-boiling behavior, including the case of one degree of freedom. The algorithm becomes more robust with increasing N-S because the possibility of finding all stationary points of the tangent-plane-distance function increases. However, the number of iterations required tends to increase with increasing N-S because the algorithm with more sampling compositions may take more iterations for merging and adding some of the sampling compositions. The general condition presented for narrow-boiling behavior is that the interplay between energy balance and phase behavior is significant. Two subsets of the condition are derived by analyzing the convex function whose gradient vectors consist of the Rachford-Rice equations; (i) the overall composition is near an edge of composition space, and (ii) the solution conditions (temperature, pressure, and overall composition) are near a critical point, including a critical endpoint. A special case of the first specific condition is the fluids with one degree of freedom. These conditions for narrow-boiling behavior are demonstrated in case studies. (C) 2016 Elsevier B.V. All rights reserved.