How can I find the eigenvalues and eigenvectors of \begin{align} Ay(p):=\int_{0}^{\infty} k^2 \cos(pk)y(k)dk \end{align} $A$ is a Hilbert-Schmidt operator.
Well actually, i came across this in solving a non-linear Hammerstein equation where the non-linearity was in the reciprocal. I am solving it numerically by first determining the Gaussian process with sample paths defined by; \begin{equation} y(k)=\sum_{i=1} ^{\infty}\frac{\alpha_i}{\sqrt{\rho_i}}u_i(k),\quad\quad y\quad\epsilon \quad L^2 \end{equation} $L^2$ is Hilbert space of functions $y$ on $[0,\Lambda]$ with the usual norm.
Also $\alpha_i$ is a sequence of random variable with mean 0 and variance 1. Here $\rho_i$ and $u_i(k)$ are eigenvalue and eigenvectors of the co-variance function defined earlier. The upper limit of the operator equation could be taken to be some cutoff $\Lambda$. Remember, the co-variance or kernel$(k^2 cos(pk))$ is non-degenerate.
It depends on the space on which you define $A$, but if we take for example $C^{\infty}\left([0,\infty[\right)$, we can choose a basis of functions and see if we can find coefficients such that a linear combination of basis functions is an eigenfunction. This is actually really similar to the finite dimensional case.
Concretely, you could use the Fourier basis $\{e^{2\pi i lx} | l \in \mathbb{Z}\}$, in which all functions of $C^{\infty}\left([0,\infty[\right)$ can be written down as $$ f(x) = \sum_{k=-\infty}^{\infty} c_k e^{2\pi i kx}. $$ Plugging this into the equation (and assuming everything converges) leads to equations for the coefficients that are solvable.