Abstract <p>A high-performance numerical method for modeling the steady-state temperature field of a heterogeneous medium is considered. A technique for determining the effective thermal conductivity of porous media is presented. The developed algorithms can be classified as collocation schemes of boundary element methods. They rely on expanding the desired solution into a series in terms of previously calculated analytical solutions to heat equations. Based on this approach, the temperature and heat flux density components at any point of the considered medium and macroscopic material parameters required for engineering applications can be computed with high accuracy and modest computer resources. The applicability and accuracy of the solution in the case of expansion functions in the form of single- and double-layer potentials are analyzed. The dependence of the effective thermal conductivity on the volumetric content of randomly distributed thermally insulated pores and spherical inclusions with thermal conductivity properties different from those of the surrounding material is studied. The algorithm is verified using well-known analytical solutions from other authors.</p>