A collocation meshless method is developed for the numerical solution of Partial Differential Equations (PDEs) on the scattered point distribution. The meshless shape functions are constructed on a group of selected nodes (stencil) arbitrarily distributed in a local support domain by means of a polynomial interpolation. This shape function formulation possesses the Kronecker delta function property, and hence many numerical treatments are as simple as those of the Finite Element Method (FEM). Nearest neighbor algorithm is used for the support domain nodes collection and a search algorithm based on the Gauss-Jordan pivot method is applied to select a suitable stencil for the construction of the shape functions and their derivatives. This search technique is subject to a monitoring procedure which selects appropriate stencil in order to keep the condition number of the resulting linear systems small. Various meshless collocation schemes for the solution of elliptic, parabolic and hyperbolic PDEs are investigated for the proposed method. Different types of PDEs are studied as test cases and all of the computational results are examined.