Non-isothermal laminar flow of a viscoelastic fluid including viscous dissipation through a square cross-section duct is analyzed. Viscoelastic stresses are described by Giesekus modele orthe Phan-Thien-Tanner model and the solvent shear stress is given by the linear Newtonian constitutive relationship. The flow through the tube is governed by the conservation equations of energy, mass, momentum associated with to one non-affine rheological model mentioned above. The mixed type of the governing system of equations (elliptic-parabolic-hyperbolic) requires coupling between discretisation methods designed for elliptic-type equations and techniques adapted to transport equations. To allow appropriate spatial discretisation of the convection terms, the system is rewritten in a quasi-linear first-order and homogeneous form without the continuity and energy equations. With the rheological models of the Giesekus type, the conformation tensor is by definition symmetrical and positive-definite, with the PTT model the hyperbolicity condition is subject to restrictions related to the rheological parameters. Based on this hyperbolicity condition, the contribution of the hyperbolic part is approximated by applying the characteristic method to extract pure advection terms which are then discretized by high ordre schemes WENO and HOUC. The algorithm thus developed makes it possible, to avoid the problems of instabilities related to the high Weissenberg number without the use of any stabilization method. Finally, a Nusselt number analysis is given as a function of inertia, elasticity, viscous dissipation, for constant solvent viscosity ratio and constant material and rheological parameters.