In this work, the stability of a two-phase, natural circulation circuit is analyzed, using a specially developed model. This thermohydraulic model results in a set of coupled, nonlinear, first-order partial differential equations, which are solved by means of the up-wind finite difference method, using combinations of explicit and implicit methods for the numerical integration of the different balance equations. An adaptive nodalization scheme is implemented, minimizing the error of the propagation of small perturbations through the discretized volumes and especially the ones having two-phase flow regime. A linearization method is implemented by means of numerical perturbations. Frequency domain calculations are carried out, allowing a rapid visualization of the stability of the linearized system. Two cases are analyzed: a test case, where the code is compared in a wide range of qualities with an analytical model, and an application case, where the model is used to analyze the stability of an integral reactor cooled by natural circulation. The CAREM prototype is taken as a reference. In both cases, the numerical diffusion and integration errors are analyzed in the stability limit prediction by means of a convergence analysis using different nodalization and numerical integration criteria.