Understanding the properties of thermal equilibrium is a fundamental problem in both classical and quantum mechanical systems. Interesting features of thermal equilibrium typically appear in systems with infinitely many degrees of freedom, presenting both intriguing physics and computational challenges. In this talk, we introduce a natural bootstrap approach to studying thermal equilibrium in classical and quantum statistical mechanics, yielding precise and rigorous numerical bounds on the observables in such systems. We apply this framework to three prototypical examples: (1) the classical Ising model on a lattice, (2) the quantum transverse field Ising model, and (3) Large N matrix quantum mechanics, which provides a window into black hole thermodynamics in the quantum gravity regime.