We propose a statistical approach to the first order phase transitions problem. For this purpose a self-consistent method for the calculation of high-order virial terms together with efficient numeric realization are developed. The method is based on the idea of differential conditions for thermodynamic consistency what makes possible the calculation of up to 30 virial terms. It is shown numerically that the obtained expansion is convergent. In the accordance with the Lee and Yang theorem, the convergency limit of virial series is found to correspond to the condensation curve. The unknown properties of high order virial terms are revealed. We show that near the critical temperature, the remainder of virial series vanishes, with high order virial terms being close to zero and the properties of a fluid being determined by the first 4–5 terms. The results are tested for two model systems; the square well and Lennard-Jones potentials. Our estimates for the phase transition curve are in good agreement with independent numeric simulation data.