Based on the nonlocal nonlinear Schrdinger equation and Poisson equation of thermal diffusion, using the slip-step Fourier algorithm and multi-grid method, we numerically investigated the propagation properties of Hermite-Guassian beams in the nonlocal thermal media. The results show that low-order Hermite-Gaussian beams can propagate stably, in contrast with the unstable propagation of high-order Hermite-Gaussian beams. The worse the stability is, the higher the order is. The effect of the boundary of the sample with different cross sections on the propagation properties of Hermite-Guassian beam is also discussed in detail. We found that propagation properties in square geometry are in agreement with those in Snyder-Mitchell model. However, in rectangular sample, the evolution of intensity distribution of Hermite-Gaussian beams differs seriously from that in the square sample.