One of the ongoing issues with fractional diffusion models is the design of an efficient high order numerical discretization. This is one of the reasons why fractional diffusion models are not yet more widely used to describe complex systems. In this paper, we derive a radial basis functions (RBF) discretization of the one-dimensional space-fractional diffusion equation. In order to remove the ill-conditioning that often impairs the convergence rate of standard RBF methods, we use the RBF-QR method [1,33]. By using this algorithm, we can analytically remove the ill-conditioning that appears when the number of nodes increases or when basis functions are made increasingly flat. The resulting RBF-QR-based method exhibits an exponential rate of convergence for infinitely smooth solutions that is comparable to the one achieved with pseudo-spectral methods. We illustrate the flexibility of the algorithm by comparing the standard RBF and RBF-QR methods for two numerical examples. Our results suggest that the global character of the RBFs makes them well-suited to fractional diffusion equations. They naturally take the global behavior of the solution into account and thus do not result in an extra computational cost when moving from a second-order to a fractional-order diffusion model. As such, they should be considered as one of the methods of choice to discretize fractional diffusion models of complex systems.