In this paper we present a new adaptive two-stage algorithm for solving elliptic partial differential equations via a radial basis function collocation method. Our adaptive meshless scheme is based at first on the use of a leave-one-out cross validation technique, and then on a residual subsampling method. Each of phases is characterized by different error indicators and refinement strategies. The combination of these computational approaches allows us to detect the areas that need to be refined, also including the chance to further add or remove adaptively any points. The resulting algorithm turns out to be flexible and effective through a good interaction between error indicators and refinement procedures. Several numerical experiments support our study by illustrating the performance of our two-stage scheme. Finally, the latter is also compared with an efficient adaptive finite element method.