  It is well known since 1960s that by exploring the tensor product structure
of the discrete Laplacian on Cartesian meshes, one can develop a simple direct
Poisson solver with an $\mathcal O(N^{\frac{d+1}d})$ complexity in d-dimension,
where N is the number of the total unknowns. The GPU acceleration of
numerically solving PDEs has been explored successfully around fifteen years
ago and become more and more popular in the past decade, driven by significant
advancement in both hardware and software technologies, especially in the
recent few years. We present in this paper a simple but extremely fast MATLAB
implementation on a modern GPU, which can be easily reproduced, for solving 3D
Poisson type equations using a spectral-element method. In particular, it costs
less than one second on a Nvidia A100 for solving a Poisson equation with one
billion degree of freedoms. We also present applications of this fast solver to
solve a linear (time-independent) Schr\"odinger equation and a nonlinear
(time-dependent) Cahn-Hilliard equation.
