  In this paper, we propose a fast and accurate numerical method based on
Fourier transform to solve Kolmogorov forward equations of symmetric scalar
L\'evy processes. The method is based on the accurate numerical formulas for
Fourier transform proposed by Ooura. These formulas are combined with
nonuniform fast Fourier transform (FFT) and fractional FFT to speed up the
numerical computations. Moreover, we propose a formula for numerical indefinite
integration on equispaced grids as a component of the method. The proposed
integration formula is based on the sinc-Gauss sampling formula, which is a
function approximation formula. This integration formula is also combined with
the FFT. Therefore, all steps of the proposed method are executed using the FFT
and its variants. The proposed method allows us to be free from some special
treatments for a non-smooth initial condition and numerical time integration.
The numerical solutions obtained by the proposed method appeared to be
exponentially convergent on the interval if the corresponding exact solutions
do not have sharp cusps. Furthermore, the real computational times are
approximately consistent with the theoretical estimates.
