  The Fokker-Planck equations for stochastic dynamical systems, with
non-Gaussian $\alpha-$stable symmetric L\'evy motions, have a nonlocal or
fractional Laplacian term. This nonlocality is the manifestation of the effect
of non-Gaussian fluctuations. Taking advantage of the Toeplitz matrix structure
of the time-space discretization, a fast and accurate numerical algorithm is
proposed to simulate the nonlocal Fokker-Planck equations, under either
absorbing or natural conditions. The scheme is shown to satisfy a discrete
maximum principle and to be convergent. It is validated against a known exact
solution and the numerical solutions obtained by using other methods. The
numerical results for two prototypical stochastic systems, the
Ornstein-Uhlenbeck system and the double-well system are shown.
