An analytical polynomial formulation of the exponential of fully symmetric matrices is proposed. This polynomial formula is derived using the Cayley - Hamilton theorem and the polynomial Euclidean division. The application field of these developments is time-dependent wavepacket calculations involving multiple potential energy surface crossing. The use of this polynomial formulation in time-dependent calculations is directly connected to the choice of the split operator, proposed by Feit and Fleck (1982 J. Comput. Phys. 47 412) as a numerical propagation scheme. We show in this paper that this analytical polynomial formula is numerically more efficient than the scheme proposed by Almeida and Metiu (1988 Chem. Phys. Lett. 146 47) where the exponential of the kinetic operator is evaluated in the momentum and diabatic representation and the exponential of the potential operator is evaluated in the coordinate and adiabatic representation. The use of an optical potential ( negative imaginary potential) to absorb the wavepacket near the edges of the grid point is quite simple to implement as it appears as a multiplicative factor in front of the polynomial expression. Our scheme is exemplified on the photodissociation of a diatomic molecule involving three potential energy curves.