The Boundary Integral Method (BIM) has recently become quite popular because of its ability to provide cheap numerical solutions to the Laplace equation. This paper describes an attempt to apply a similar approach to the (time-dependent) heat equation in two space variables. 1. Introduction. In its simplest form, the BIM uses the prescribed initial and boundary data, together with the (known) fundamental solution of a given differen- tial equation defined on some domain Q, to construct a second integral equation which is itself to be solved on the boundary of Ui. In our case, the integral equation solution is a layered thermal potential, i.e., it represents a continuous distribution of sources and sinks along the boundary of the two-dimensional domain U. One then obtains the solution to the given problem by integrating the product of the layered thermal potential and a known kernel over the boundary, hence the terminology boundary integral. The advantage of such an approach is that the heart of the computation, viz. the solution of the integral equation, is performed on the boundary, thus reducing the problem from two space dimensions to one. Furthermore, the resulting domain in one space and time will be rectangular, a computational convenience, even though the given domain i2 may have been quite irregular (see Example 5 in the Appendix). Of course, as a result of these simplifications, we may expect substantial savings in computer time and storage. The work outlined below is based on the use of single layered thermal potentials, and requires that the domain have a smooth boundary and, with some restrictions, either a Neumann or mixed boundary condition. The equation itself must be homogeneous, but we do allow inhomogeneous initial data. Though not shown here, it seems to be well within the capabilities of this approach to handle boundaries involving an arbitrary mixture of piecewise C2 curves, with Dirichlet, Neumann, and/or mixed boundary data. Of course, the tradeoff is that to evaluate the solution, we must do a double integral for each point at which we want to know the solution. If needed at a large number of points, the cost of generating the solution dominates, and the method becomes impractical. In many applications however, in particular in semiconductor