Computations of incompressible flows with all velocity boundary conditions require solution of a Poisson equation for pressure or pressure-corrections with all Neumann boundary conditions. Discretization of such a Poisson equation results in a rank-deficient ill-conditioned matrix of coefficients. When a non-conservative discretization method such as finite difference, finite element, or spectral scheme is used, the ill-conditioned matrix also generates an inconsistency which makes the residuals in the iterative solution to saturate at a threshold level that depends on the spatial resolution and the order of the discretization scheme. In this paper, we examine this inconsistency for a high-order meshless discretization scheme suitable for solving the equations on a complex domain. The high order meshless method uses polyharmonic spline radial basis functions (PHS-RBF) with appended polynomials to interpolate scattered data and constructs the discrete equations by collocation.