We present a theoretical analysis of a fictitious domain formulation of the Newtonian cooling problem, motivated by applications in topology optimization. The method reformulates the classical heat conduction model with Robin-type boundary conditions on a fixed computational domain using a so-called weak material approximation. In this setting, the conductivity equals one in the solid subdomain Omega s and a small positive parameter & varepsilon; in the surrounding fictitious region. We derive a priori error estimates that quantify the consistency error between the extended and original formulations and prove that the solution restricted to Omega s converges to the true solution with an O(& varepsilon;) error in the H1(Omega s) norm. We further provide & varepsilon; -dependent finite element (FE) error estimates and show that for a mesh with characteristic size h, the condition number of the FE systems scales as O(& varepsilon;-1h-2) . To address the & varepsilon; -induced ill-conditioning, we introduce a simple yet effective diagonal preconditioning strategy that removes the & varepsilon; -dependency from the condition number. Numerical experiments confirm the theoretical convergence rate and demonstrate the effectiveness of the method. Thus, this work provides a theoretical foundation for utilizing weak material approximations for boundary-effect-dominated problems, thereby extending existing analyses to cases with Robin-type boundary conditions.