A new finite element scheme for the numerical simulation of 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 spatial coordinate system attached to the particles is discretized by ten-node quadratic tetrahedral elements using Cartesian coordinates and the pressure by linear interpolation inside these elements. The spatial discretization of the governing equations follows the mixed Galerkin finite element method. The time integral is discretized by a quadratic interpolation in time. The convergence of the method in time and space was demonstrated on the free surface problem of a filament stretched between two plates, considering the axisymmetric case as well as the growth of non-axisymmetric disturbances on the free surface. The scheme converges with respect to discretization of the spatial and time dimensions, both with third order.