A numerical method for the study of the gravothermal instability in star clusters