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 \(\epsilon \) 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(\epsilon )\) error in the \(H^1(\Omega _s)\) norm. We further provide \(\epsilon \) -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(\epsilon ^{-1}h^{-2})\) . To address the \(\epsilon \) -induced ill-conditioning, we introduce a simple yet effective diagonal preconditioning strategy that removes the \(\epsilon \) -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.