Krylov subspace methods are a popular class of iterative algorithms for solving Ax=b via matrix-vector products, with their convergence rates often depending on the conditioning and normality of A. In this talk we develop analogues of the conjugate gradient method, MINRES, and GMRES for the spectral solution of Lu = f via operator-function products, where L is an unbounded elliptic differential operator. Our Krylov methods employ operator preconditioning and have good convergence properties. This setting allows us to automatically preserve the continuous properties of the operator regardless of the basis chosen, and leads to competitive solvers compared to classical spectral discretizations.