A new numerical approach for the simulation of non-isothermal three-dimensional time-dependent flow of viscoelastic fluids is presented. The viscoelastic fluids are of the K-BKZ integral type and the method is based on a Lagrangian kinematics description of the fluid flow. The K-BKZ fluid is assumed to be a thermorheological simple material using the extended Morland and Lee hypothesis by Crochet and Naghdi [M.J. Crochet, RM. Naghdi, A class of non-isothermal viscoelastic fluids, International journal of Engineering Science 10 (1972) 755-800], where the real time in the K-BKZ constitutive equation is replaced with a temperature dependent pseudo time. The spatial coordinate system attached to the particles is discretized by 10-node quadratic tetrahedral elements using Cartesian coordinates. The temperature and the pressure are discretized by 10-node quadratic and linear interpolation, respectively, in the tetrahedral particle elements. The spatial discretization of the governing equations follows a mixed Galerkin finite element method. This type of scheme ensures third order accuracy with respect to the discretization of spatial dimension. The temperature equation is solved in time utilizing an implicit variable step backwards differencing (BDF2) scheme, obtaining second order convergence of the temperature in time. A quadratic interpolation in time is applied to approximate the time integral in the K-BKZ equation. This type of scheme ensures third order accuracy with respect to the discretization of the time. The third order convergence in time and space is demonstrated on time-dependent non-isothermal elongational flow as a test case. (C) 2009 Published by Elsevier B.V.