The thermal boundary layer flow is a canonical flow with characteristics that are present in most natural and industrial convection flows. An approximate self-similar solution is proposed for the first time for the thermal boundary layer of steady laminar flow of viscoelastic fluids, described by the finitely extensible nonlinear elastic constitutive equation with Peterlin's closure (FENE-P model). This semi-analytical thermal solution is obtained by performing an order of magnitude analysis and ensuing simplifications of the governing equations by assuming that the fluid properties are independent of temperature therefore decoupling the flow governing equations from the energy equation. The effects of viscoelasticity quantified with the Weissenberg number based on the streamwise coordinate (x) (Wix) up to Wix=1 and viscous dissipation (results are presented for Brinkman numbers between -40 and +40) on thermal boundary layer characteristics are investigated comprehensively for both constant wall temperature and constant wall heat flux. At low elasticity levels (Wix<0.01) the solution exhibits a global self-similar behavior in which flow and thermal quantities collapse on the corresponding Newtonian curves, and the polymer characteristics show a unique behavior if adequately normalized. However, by increasing flow elasticity the unique self-similar behavior of the approximate solution is lost, with the elasticity dependent results exhibiting local variations. In addition, the effects of elasticity are intensified by viscous dissipation. For the present study cases, it is observed that elasticity may change Nusselt numbers by more than 8%, and the thermal boundary layer thickens by up to 10%.