An iterative boundary integral numerical method for solving the steady conduction of heat is developed. The method is general for two- and three-dimensional regions with arbitrary boundary shapes. The development is generalized to include the first, second, and third kind of boundary conditions and also radiative boundary and temperature-space dependent convective coefficient cases. With Kirchhoff’s transformation, cases of temperature-dependent thermal conductivity with general boundary conditions are also accounted for by the present method. A variety of problems are analyzed with this method and their solutions are compared to those obtained analytically. A comparison between the present method and the finite difference predictions is also investigated for a case of mixed temperature and convective boundary conditions. Moreover, two-dimensional regions with three kinds of boundary conditions and irregular-shaped boundaries are used to illustrate the versatility of the technique as a computational procedure.