The solution of the Helmholtz equation is a fundamental step in frequency domain seismic imaging. This paper deals with a numerical study of solutions for 2D Helmholtz equation using a Gaussian radial basis function-generated finite difference scheme (RBFFD). We analyze the behavior of the local truncation error in approximating partial derivatives of the 2D Helmholtz equation solutions when the shape parameter of RBF varies. For discretization, we performed, by means of a classical numerical dispersion analysis with plane waves, a minimization of the error function to obtain local and adaptive near optimal shape parameters according to the local wavelength of the required solution. In particular, the method is applied to obtain a simple and accurate solver by using stencils which seven nodes on hexagonal regular grids, wich mitigate pollution-effects. We validated numerically that the stability and isotropy are enhanced with respect to Cartesian grids. Our method is tested with standard case studies and velocity models, showing similar or better accuracy than finite difference and finite element methods. This is an efficient way for interacting with inverse and imaging problems such as Full Wave Inversion.