Makino (1986) performed 100-body simulations of an equal-mass cluster, using both softened and unsoftened potentials. They used a softened potential in order to reduce the energy production from binaries. They found oscillatory behavior for both cases, but with larger amplitude for softened potential. They argued that the fact that the amplitude of oscillation is larger for reduced energy production implies that the observed oscillation is driven by the gravothermal instability and not purely by the heating by binary. McMillan and Lightman (1984, [@McMillan1986a]) developed hybrid calculation code in which the central region of a cluster is treated as an $N$-body system while the outer region is treated in FP formalism. They observed oscillatory behavior but concluded that it was not gravothermal. Inagaki (1986) performed 1000- and 3000-body simulations but found no sign of oscillation. He used standard Salpeter mass function as the distribution of the masses of particles. Makino (1989) performed simulation of a 3000-body equal-mass system and saw some signs of oscillation, including the temperature inversion. Spurzem and Aarseth (1996) performed 10,000 body simulation and saw oscillatory behavior, but they did not see any clear evidence of the gravothermal oscillation. $N$-body simulations in the last decade can be summarized 