A boundary integral numerical method of solving the steady conduction of heat in composite media is developed. The method is suitable for solving two and three dimensional regions with arbitrary boundary shapes. The present method is general for cases with first, second and third boundary conditions. A variety of problems are analyzed using the method and presented for illustration of its use and capabilities. A comparison between the present method and the finite difference predictions is also made. Finally, a composite region with the three kinds of boundary conditions present and irregular shaped boundaries is used to demonstrate the versatility of the technique. The method is proved to be efficient and accurate as a computational procedure for solving the general problem of heat conduction in composite media.