Multiple steady-state solutions have been found theoretically for a two-phase natural circulation loop. This bifurcation phenomenon is attributed to the nonmonotonic behavior of the buoyancy and friction forces, and the discontinuities due to transitions between two-phase flow patterns. The paper outlines a general and consistent modeling method to describe two-phase natural circulation at steady state. The solution of the one-dimensional governing conservation equations is performed by a combination of analytical and numerical stages. The method is implicit, to treat the highly nonlinear problem by multinested iteration loops. Results are presented and discussed for the multiple solutions of the flow rate, temperature, quality and void distributions and two-phase elevations. Pressure effects are also studied. A comparison with available data shows that the theoretical model simulates the important characteristics of two-phase thermosyphons.