We present an efficient integral equation approach to solve the forced heat equation, u t ( x ) − Δ u ( x ) = F ( x , u , t ) , in a two-dimensional, multiply-connected domain, with Dirichlet boundary conditions. Instead of using an integral equation formulation based on the heat kernel, we discretize in time, first. This approach, known as Rothe’s method, leads to a non-homogeneous modified Helmholtz equation that is solved at each time step. We formulate the solution to this equation as a volume potential plus a double layer potential, and both of these potentials are calculated with available tools accelerated by the fast multipole method. For a total of N points in the discretization of the boundary and the domain, the total computational cost per time step is O ( N ) . We demonstrate our approach on the heat equation and the Allen–Cahn equation.