Radial basis function (RBF) approximation is an extremely powerful tool for representing smooth functions in non-trivial geometries, since the method is meshfree and can be spectrally accurate. A perceived practical obstacle is that the interpolation matrix becomes increasingly ill-conditioned as the RBF shape parameter becomes small, corresponding to flat RBFs. Two stable approaches that overcome this problem exist, the Contour-Padé method and the RBF-QR method. However, the former is limited to small node sets and the latter has until now only been formulated for the surface of the sphere. This paper contains an RBF-QR formulation for planar two-dimensional problems. The algorithm is perfectly stable for arbitrarily small shape parameters and can be used for up to a thousand node points in double precision and for several thousand node points in quad precision. A sample MATLAB code is provided.