This study deals with the case of multiple internal heat source inversion problem of steady-state based on nondestructive infrared detection. We construct homogeneous and heterogeneous steady heat conduction models of different shapes. Neither the number of heat sources, nor their locations, nor their sizes nor their intensities are known. We use the finite element method (FEM) based on numerical algorithm to analyze the two-dimensional model discretely. The internal heat conduction process of model is analyzed. The resultant temperature field can be decomposed into the temperature field caused by the ambient temperature and those given by internal heat sources. We simplify the finite element matrix equation according to decomposition process above. Finally, the problem boils down to solving the highly underdetermined matrix equation of Ax = b. The unknown x item corresponds to internal thermal heat source field. The piecewise polynomial spectral truncated singular value decomposition (PPTSVD) is applied for the first time to the inverse heat source problem. Its regular operator matrix is changed from the original more order differential operator matrix to regional node weighted matrix. After replacement this solution improves the effect of heat source field tending to the boundary. Results of the solution confirms a real heat source field when there are less heat sources or different heat sources are far from each other. But there also exits a serious superimposed effect between neighboring heat sources. We improve the algorithm to study this problem through using the iterative elimination process which complies with the idea of spreading heat source field and then gathering. The iteration tolerance and number of times belonging to one single PPTSVD solving process are reduced. Through iterating the multiple PPTSVD solving process and reconstructing matrixes A and b in each iteration, we obtain the scatter heat source field distribution surrounding real field. Finally, this scattered distribution solution is gathered again. According to the heat source parameter calculated by algorithm, the temperature field of whole model can be reconstructed using FEM. Comsol numerical simulation and real physical experiments are performed to verify the validity and accuracy of the algorithm in different heat conduction models. The results demonstrate that the algorithm can access each parameter of multiple heat sources. Even in the heterogeneous model, it can still obtain accurate results and reconstruct the temperature field of the two-dimensional model. The algorithm can be applied to non-destructive material detection and human infrared medical imaging fields.