A major obstacle to the application of the standard Radial Basis
Function-generated Finite Difference (RBF-FD) meshless method is constituted by
its inability to accurately and consistently solve boundary value problems
involving Neumann boundary conditions (BCs). This is also due to
ill-conditioning issues affecting the interpolation matrix when boundary
derivatives are imposed in strong form. In this paper these ill-conditioning
issues and subsequent instabilities affecting the application of the RBF-FD
method in presence of Neumann BCs are analyzed both theoretically and
numerically. The theoretical motivations for the onset of such issues are
derived by highlighting the dependence of the determinant of the local
interpolation matrix upon the boundary normals. Qualitative investigations are
also carried out numerically by studying a reference stencil and looking for
correlations between its geometry and the properties of the associated
interpolation matrix. Based on the previous analyses, two approaches are
derived to overcome the initial problem. The corresponding stabilization
properties are finally assessed by succesfully applying such approaches to the
stabilization of the Helmholtz-Hodge decomposition