  The quadratic numerical range $W^2(A)$ is a subset of the standard numerical
range of a linear operator which still contains its spectrum. It arises
naturally in operators which have a $2 \times 2$ block structure, and it
consists of at most two connected components, none of which necessarily convex.
The quadratic numerical range can thus reveal spectral gaps, and it can in
particular indicate that the spectrum of an operator is bounded away from $0$.
  We exploit this property in the finite-dimensional setting to derive Krylov
subspace type methods to solve the system $Ax = b$, in which the iterates arise
as solutions of low-dimensional models of the operator whose quadratic
numerical ranges is contained in $W^2(A)$. This implies that the iterates are
always well-defined and that, as opposed to standard FOM, large variations in
the approximation quality of consecutive iterates are avoided, although $0$
lies within the convex hull of the spectrum. We also consider GMRES variants
which are obtained in a similar spirit. We derive theoretical results on basic
properties of these methods, review methods on how to compute the required
bases in a stable manner and present results of several numerical experiments
illustrating improvements over standard FOM and GMRES.
