A popular approach for modeling and inference in spatial statistics is to represent Gaussian random fields as solutions to stochastic partial differential equations (SPDEs) of the form (Formula presented.), where (Formula presented.) is Gaussian white noise, L is a second-order differential operator, and (Formula presented.) is a parameter that determines the smoothness of u. However, this approach has been limited to the case (Formula presented.), which excludes several important models and makes it necessary to keep β fixed during inference. We propose a new method, the rational SPDE approach, which in spatial dimension (Formula presented.) is applicable for any (Formula presented.), and thus remedies the mentioned limitation. The presented scheme combines a finite element discretization with a rational approximation of the function (Formula presented.) to approximate u. For the resulting approximation, an explicit rate of convergence to u in mean-square sense is derived. Furthermore, we show that our method has the same computational benefits as in the restricted case (Formula presented.). Several numerical experiments and a statistical application are used to illustrate the accuracy of the method, and to show that it facilitates likelihood-based inference for all model parameters including β. Supplementary materials for this article are available online.