The paper presents results on the approximation of functions which solve an elliptic differential equation by operator adapted systems of functions. Compared with standard polynomials, these operator adapted systems have superior local approximation properties. First, the case of Laplace's equation and harmonic polynomials as operator adapted functions is analyzed and rates of convergence in a Sobolev space setting are given for the approximation with harmonic polynomials. Special attention is paid to the approximation of singular functions that arise typically in corners. These results for harmonic polynomials are extended to general elliptic equations with analytic coefficients by means of the theory of Bergman and Vekua; the approximation results for Laplace's equation hold true verbatim, if harmonic polynomials are replaced with generalized harmonic polynomials. The Partition of Unity Method is used in a numerical example to construct an operator adapted spectral method for Laplace's equation that is based on approximating with harmonic polynomials locally.