A novel temporal finite element method is presented to solve static viscoelastic problems, which is more flexible and appropriate to describe temporal variation of displacement than conventional finite difference or numerical integral methods. By using virtual work principle, a spatial finite element based governing equation is derived in terms of displacement and its derivatives. Then two temporal finite element models are developed using Gurtin variational principle and weighted residual technique. A kind of hybrid shape functions with polynomial and trigonometric basis is stressed to give more flexible descriptions of time varying variables. A criterion of stability analysis is derived, which can numerically be conducted when the constitution of shape functions and step size are prescribed. A recursive algorithm toward end is developed by which the temporal FE analysis can be conducted via a matrix power product with the initial condition, instead of step marching. The proposed approach is available for the viscoelastic model described by linear differential equations with order h ≤ 2, and can conveniently be combined with well-developed numerical algorithms to tackle with boundary value problems, such as FEM, SBFEM etc. Various numerical examples, including those with static/harmonic loads, creep, stress singularity and heterogeneous structures, etc. are provided to illustrate the efficiency of the proposed approaches, and impacts of the temporal FE model, constitution of shape functions, and step size, etc. are taken into account. Satisfactory results are achieved in comparison with analytical or ABAQUS-based solutions.