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)

<b>Edit:</b>A number of commentators have described the natural technique:  discretize space, turn the covariance into a matrix, easily solve for the eigenvectors, then reconstruct the functions $f_k(t)$ by interpolation.  

Our motivation is to generate a Gaussian random field $\xi$ using the Karhunen-Loève transformation:  let $Z_1, \cdots, Z_K$ be independent Gaussian RVs, and define the random field $$\xi(t) = \sum_{k=1}^K Z_k f_k(t).$$  Since the covariance is so smooth, the eigenfunctions will be smooth, and the Gaussian field will be smooth.  I'm sure that a competent numerical analyst could use the discretization technique to reconstruct the eigenfunctions up to arbitrary derivatives.  Unfortunately, I am an incompetent numerical analyst, and I would surely introduce systematic errors along the way.  Yet the problem is so common that I am still hoping somebody could point us to a simple reference which does this well.

To make the question more mathematical precise:  the numerical approximations $f_k$ should be within an arbitrarily small $C^3$-error of the actual eigenfunctions.