Abstract. This article describes the derivation and implementation of a numerical method to solve constant-coefficient, parabolic partial differential equations in two space dimensions on rectangular domains. The method is based on a formula for the Green's function for the problem obtained via reflections at the boundary of the domain from the corresponding formula for the fundamental solution in the whole plane. It is inspired by a related method for variable coefficients equations in the whole space introduced by Constantinescu, Costanzino, Mazzucato, and Nistor in J. Math. Phys, 51 103502 (2010). The benchmark case of the two-dimensional heat equation is considered. We compare the Green's function method with a finite-difference scheme, more precisely, an alternating direction implicit (ADI) method due to Peaceman and Rachford. Our method yields better rates of convergence to the exact solution.