Non-isothermal laminar flow of a viscoelastic fluid through a square cross-section duct is analyzed. Viscoelastic stresses are described by the Phan-Thien -Tanner model and the solvent shear stress is given by the linear Newtonian constitutive relationship. The solution of the set of governing equations spawns coupling between equations of elliptic-hyperbolic type. Our numerical approach is based on the finite-differences method. To treat the hyperbolic part, the system of equations are rewritten in a quasilinear form. The resulting pure advection terms are discretized using high-order upwind schemes when the hyperbolicity condition is satisfied. The incompressibility condition is obtained by the semi-implicit projection method. Finally we investigate the evolution of velocity, shear stress, viscosity and heat transfer over a wide range of Weissenberg numbers.