  We present a practical algorithm to approximate the exponential of
skew-Hermitian matrices based on an efficient computation of Chebyshev
polynomials of matrices and the corresponding error analysis. It is based on
Chebyshev polynomials of degrees 2, 4, 8, 12 and 18 which are computed with
only 1, 2, 3, 4 and 5 matrix-matrix products, respectively. For problems of the
form $\exp(-iA)$, with $A$ a real and symmetric matrix, an improved version is
presented that computes the sine and cosine of $A$ with a reduced computational
cost.
  The theoretical analysis, supported by numerical experiments, indicates that
the new methods are more efficient than %diagonalising the matrix when its norm
is not very large, and also than other schemes based on rational Pad\'e
approximants and Taylor polynomials for all tolerances and time interval
lengths. The new procedure is particularly recommended to be used in
conjunction with exponential integrators for the numerical time integration of
the Schr\"odinger equation.
