A numerical integration method is presented for the treatment of transient heat conduction problems. A Cartesian formulation is developed that is suitable for the treatment of irregular domains under general boundary conditions. The qualities of the scheme are demonstrated, in terms of both accuracy and computational efficiency, by comparison with analytical and numerical solutions. Results for the basic two-dimensional annular geometry show that the method has nearly second-order accuracy in space and time, at least in simple cases. Finally, a complex multiconnected domain is considered, to test the method performance under more severe conditions, including the presence of multiple length scales. The numerical experiment demonstrates that the numerical scheme is efficient, stable, and convergent.