A numerical mathematical model was developed to determine the behavior of thermal distribution in a convergent-divergent nozzle used in aerospace applications. This development was carried out by using the two-dimensional heat equation for describing the phenomenon, which was resolved by means of the finite difference method. In addition, it was considered insulated and flow boundary conditions, a structured mesh with, and a nodal density of 169. The methodology was implemented in MATLAB through an algorithm which allows variations in the geometry, nodal density, and inlet conditions. This methodology was validated using computational fluid dynamics (CFD). The validation demonstrated that the proposed methodology reduces the computational cost by having an execution time 1555 times less than CFD simulation time. The relative error between the methodology and the CFD results was 2,9 %. The proposed methodology would allow the preliminary design process of this type of nozzles by analyzing thermal distribution effectively and accurately.