The aim of this paper is to provide tools and results for the analysis of the linear systems arising from radial basis function (RBF) approximations of partial differential equations (PDEs), see e.g., [1,9]. Informally, a radial function phi(x) :R-n -> R is a function of the Euclidean norm parallel to x parallel to of x, i.e., phi(x) = eta(parallel to x parallel to), for n(t): R -> R. Examples are functions of the following form root t(2) + c(2), multiquadric (MQ), 1 root t(2) + c(2), inverse multiquadric (IMQ), e-(t2/c2), Gaussian. In this context c is the shape parameter, whose value plays a role in modeling problems with various specific features. At least numerically, it is evident that the precision of the approximation procedures based on RBFs is very high. In fact, if h denotes the maximal step size, then the approximation error behaves like O(lambda(c/h)) for the MQ and like O(lambda(root c/h)) for the IMQ and Gaussian, where lambda is a positive parameter, strictly less than one, and independent of h. The price that has to be paid concerns the increasing ill -conditioning of the related linear systems in which a growth of the order of e(theta c/h) is observed at least for large values of c/h, with theta being a positive constant independent of h and also of the shape parameter c. We are interested in the spectral behavior of the resulting matrices, and especially in the extremal behavior (conditioning) and in the global distribution results: such a study is crucial for designing fast and accurate solution methods. A first important step in understanding the spectral behavior of the considered matrices was done in [3], where the remarkable link with Toeplitz sequences generated by a symbol was exploited. Here we give a more precise analysis than in [3], by showing that for some choices of RBF, e.g., IMQ, Gaussian, and for some values of the parameter c/ h, the conditioning is not