  We propose a unified meshless method to solve classical and fractional PDE
problems with $(-\Delta)^{\frac{\alpha}{2}}$ for $\alpha \in (0, 2]$. The
classical ($\alpha = 2$) and fractional ($\alpha < 2$) Laplacians, one local
and the other nonlocal, have distinct properties. Therefore, their numerical
methods and computer implementations are usually incompatible. We notice that
for any $\alpha \ge 0$, the Laplacian $(-\Delta)^{\frac{\alpha}{2}}$ of
generalized inverse multiquadric (GIMQ) functions can be analytically written
by the Gauss hypergeometric function, and thus propose a GIMQ-based method. Our
method unifies the discretization of classical and fractional Laplacians and
also bypasses numerical approximation to the hypersingular integral of
fractional Laplacian. These two merits distinguish our method from other
existing methods for the fractional Laplacian. Extensive numerical experiments
are carried out to test the performance of our method. Compared to other
methods, our method can achieve high accuracy with fewer number of unknowns,
which effectively reduces the storage and computational requirements in
simulations of fractional PDEs. Moreover, the meshfree nature makes it free of
geometric constraints and enables simple implementation for any dimension $d
\ge 1$. Additionally, two approaches of selecting shape parameters, including
condition number-indicated method and random-perturbed method, are studied to
avoid the ill-conditioning issues when large number of points.
