scarce. In this regime, Gibbs ensemble Monte Carlo (GEMC) cannot reproduce binodal curves [6, 32]. The phase diagram at large κ was studied by Monte Carlo supplemented by thermodynamic integration (MC-TI) in Ref. [7], but the results do not include the critical density. In our previous works we reported the coexistence properties of HCAY fluid for long range tails. These properties were compared to those reported in the literature, and some differences were found among them [33–35]. In a recent work, Singh [36] has reported the coexistence curves of short-range attractive Yukawa (SR-HCAY) fluid with κ $^{= 8}$−10, using Grand-canonical transition-matrix Monte Carlo (GC-TMMC) with the histogram reweight- ing method. However, this method proves difficult to use in order to calculate coexistence properties at low temperature and high liquid densities [36]. So, NVT-MC simulations ap- pear to be a more efficient method to calculate such properties for very short-range systems, where these conditions are met. In view of the above considerations, the main goals of this paper are to demonstrate that there is a linear dependence between the critical density and the reciprocal of the critical temperature even for short interaction ranges, where theoretical approaches show different tendencies, and to report a systematic study of the liquid-vapor phase diagrams of SR- HCAY fluid with κ = 9 , 10 , 12 and 15, using canonical ensemble Monte Carlo simulation, where most theoretical approaches fail and other simulation techniques cannot reproduce the binodal curves. Finally, once again we have confirmed that the coexistence curves of SR- HCAY model follow the law of corresponding states, i.e., curves corresponding to different κ fall on top of each other within a high degree of accuracy, provided the density ρ and temperature T are rescaled by their critical values ρ$_{c}$ , T$_{c}$ . Besides, we have also reported the binodal curves and critical points predicted by SCOZA for the above values of κ . Although the application of SCOZA to narrow Yukawa potentials has been considered before [1, 23], an extensive comparison with simulation data for the phase diagram in this regime has not been possible so far because of the aforementioned paucity of simulations at large κ , and the present study represents a good opportunity to perform it. Results were obtained both by the standard version of SCOZA considered in Refs. [1, 23, 28], and by a modified version developed in [29], where consistency with the virial route, which is disregarded in the original implementation, is partially taken into account. The comparison shows that some discrepancies between SCOZA and simulations are definitely present: specifically, the standard formulation of SCOZA overestimates both 3