The numerical solution of partial differential equations by discretization techniques is ubiquitous in computational physics. In this work we benchmark this approach in the quantum realm by solving the heat equation for a square plate subject to fixed temperatures at the edges and random heat sources and sinks within the domain. The hybrid classical-quantum approach consists in the solution on a quantum computer of the coupled linear system of equations that result from the discretization step. Owing to the limitations in the number of qubits and their connectivity, we use the Gauss-Seidel method to divide the full system of linear equations into subsystems, which are solved iteratively in block fashion. Each of the linear subsystems were solved using 2000Q and Advantage quantum computers developed by D-Wave Systems Inc. By comparing classical numerical and quantum solutions, we observe that the errors and chain break fraction are, on average, greater on the 2000Q system. Unlike the classical Gauss-Seidel method, the errors of the quantum solutions level off after a few iterations of our algorithm. This is partly a result of the span of the real number line available from the mapping of the chosen size of the set of qubit states. We verified this by using techniques to progressively shrink the range mapped by the set of qubit states at each iteration (increasing floating-point accuracy). As a result, no leveling off is observed. However, an increase in qubits does not translate to an overall lower error. This is believed to be indicative of the increasing length of chains required for the mapping to real numbers and the ensuing limitations of hardware.