The Monte Carlo entropic sampling method previously presented in [I. Shteto, J. Linares, F. Varret, Phys. Rev. E 56 (1997) 5128] is adapted here to an Ising-like system with short- and long-range interactions. Such model is suited to spin crossover solids [J. Linares, H. Spiering, F. Varret, Eur. J. Phys. B 10 (1999) 271; K. Boukheddaden, J. Linares, H. Spiering, F. Varret, Eur. Phys. J. B 15 (2000) 317] where the long interaction is due to elastic coupling mediated by the lattice, while the short-range interaction originates from the bonding between the spin crossover units [J. Linares, H. Spiering, F. Varret, Eur. J. Phys. B 10 (1999) 271]. Taking into account the different degeneracies gHS for high-spin (HS) and gLS for low-spin (LS) states, the Ising Hamiltonian associated with fictitious spins is written: H = -h ∑ σi - J ∑ σiσj with h = -Δ/2+kBT ln(gHS/gLS/2+G〈σ〉 where J and G are the short- and long-range interactions, respectively, and A the energy gap of ligand field such that the LS state is the ground state. The numerical method has been tested successfully by comparison to the exact solution for a 1D system: [Fe(Htrz)2(trz)](BF4)2 [J. Linares, H. Spiering, F. Varret, Eur. J. Phys. B 10 (1999) 271; J. Krober, J.P. Audière, R. Claude, O. Kahn, J. Hassnoot, F. Grolière, C. Jay, A. Bousseksou, J. Linares, F. Varret, A. Gonthier-Vassal, Chem. Mater. 6 (1994) 1404]. We describe here the results obtained for 2D systems, and show that the squareness of the thermal hystersis loop, associated with the spin-transition, can be correlated to the strength of short-range interactions. © 2003 Elsevier Science Ltd. All rights reserved.