  The availability of efficient Krylov subspace solvers play a vital role for
the solution of a variety of numerical problems in computational science. Here
we consider lattice field theory. We present a new general numerical method to
compute many Green's functions for complex non-singular matrices within one
iteration process. Our procedure applies to matrices of structure $A=D-m$, with
$m$ proportional to the unit matrix, and can be integrated within any Krylov
subspace solver. We can compute the derivatives $x^{(n)}$ of the solution
vector $x$ with respect to the parameter $m$ and construct the Taylor expansion
of $x$ around $m$. We demonstrate the advantages of our method using a minimal
residual solver. Here the procedure requires $1$ intermediate vector for each
Green's function to compute. As real life example, we determine a mass
trajectory of the Wilson fermion matrix for lattice QCD. Here we find that we
can obtain Green's functions at all masses $\geq m$ at the price of one
inversion at mass $m$.
