Let $c(r)$ be a nice, continuous function with compact support.  For example, $c(r) = \tfrac 1 5 (1-r)^{11} \big( 5 + 55r + 239 r^2 + 429 r^3 \big)$ for $r \in [0,1]$, and $c(r) = 0$ otherwise.  

On a rectangular domain $D$ in $\mathbb R^2$, we define $c(x,y) = c(|x-y|)$, which we interpret as the kernel of an integral operator.  Consider the integral equation $$\iint_D c(x,y) f(y) \  \mathrm{d}y = \lambda f(x).$$  By <a href="http://en.wikipedia.org/wiki/Mercer%27s_theorem">Mercer's theorem</a>, we can find a sequence of non-negative eigenvalues $\lambda_k$ and eigenfunctions $f_k(x)$ such that:

* The sequence $f_k$ is an orthonormal basis for $L^2(D)$,

* The eigenfunctions $f_k$ corresponding to non-negative eigenvalues are continuous on $D$, and

* We have the representation $c(x,y) = \sum \lambda_k f_k(x) f_k(y)$.

<b>Question:</b>  How can we numerically approximate the first $K$ eigenfunctions $f_1, \cdots, f_K$ and keep the total error within a desired error?  (Polynomial or Fourier approximations would suffice)