An (N+2)-field model proposed in our previous study was modified so as to be applicable to compressible two-phase flows by taking into account the dependence of phase densities on temperature and pressure. Then, a numerical method for solving the modified (N+2)-field model was developed. The main feature of the proposed method is that the numerical stability is not determined by the Courant-Friedrichs-Levy condition but by the Courant condition, and therefore, stable numerical solutions can be obtained with a time step &UDelta; t much larger than the one used in a standard numerical method for solving compressible two-phase flows. The verification of the proposed method was carried out through three sample calculations, i.e., simulations of an adiabatic air-water bubbly flow in a vertical duct of 10m in length, an air-water bubbly flow in a uniformly heated duct of 10m in length and a bubbly flow in a large bubble column. As a result, it was confirmed that (a) the proposed method accurately predicts the effects of pressure and temperature on phase densities and volumetric fluxes and (b) a simulation of a compressible bubbly flow in a practical bubble column can be performed with &UDelta; t satisfying the Courant condition.