Abstract: This paper proposes an accurate and computationally efficient solver of the heat equation ( δ t + Δ ) F ( · , t ) = 0 , F ( · , 0 ) = f , on a volumetric domain, through the (r,r)-degree Padé-Chebyshev rational approximation of the exponential representation F ( · , t ) = exp ( − t Δ ) f of the solution. To this end, the heat diffusion problem is converted to a set of r differential equations, which involve only the Laplace–Beltrami operator, and whose solution converges to F ( · , t ) , as r → + ∞ . The discrete heat equation is equivalent to r sparse, symmetric linear systems and is independent of the volume discretization as a tetrahedral mesh or a regular grid, the evaluation of the Laplacian spectrum, and the selection of a subset of eigenpairs. Our approach has a super-linear computational cost, is free of user-defined parameters, and has an approximation accuracy lower than 10−r. Finally, we propose a simple criterion to select the time value that provides the best compromise between approximation accuracy and smoothness of the solution.