One of the methods for calculating time propagators in quantum mechanics uses an expansion of e(-iHt/h) in a sum of orthogonal polynomial. Equations involving Chebychev, Legendre, Laguerre, and Hermite polynomials have been used so far. We propose a new formula, in which the propagator is expressed as a sum in which each term is a Gegenbauer polynomial multiplied with a Bessel function. The equations used in previous work can be obtained from ours by giving specific values to a parameter. The expression allows analytic continuation from imaginary to real time, transforming thus results obtained by evaluating thermal averages into results pertaining to the time evolution of the system. Starting from the expression for the time propagator we derive equations for the Green's function and the density of states. To perform computations one needs to calculate how the polynomial in the Hamiltonian operator acts on a wave function. The high order polynomials can be obtained from the lower ordered ones through a three term recursion relation; this saves storage and computer time. As a numerical test, we have computed the bound state spectrum of the Morse oscillator and the transmission coefficient for tunneling through an Eckart barrier. We have also studied the evolution of a Gaussian wave packet in a Morse potential well.