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 Mercer's theorem, 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)$.

Question: 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)