In this paper, we employ linear combinations of n heat kernels to approximate solutions to the heat equation. We show that such approximations are of order O(t(1/2p - 2n+1/2)) in L(p)-norm, 1 <= p <= infinity, as t -> infinity. For positive solutions of the heat equation such approximations are achieved using the theory of truncated moment problems. For general sign-changing solutions these type of approximations are obtained by simply adding an auxiliary heat kernel. Furthermore, inspired by numerical computations, we conjecture that such approximations converge geometrically as n -> infinity for any fixed t > 0.