  We consider the problem of approximating the solution to $A(\mu) x(\mu) = b$
for many different values of the parameter $\mu$. Here we assume $A(\mu)$ is
large, sparse, and nonsingular with a nonlinear dependence on $\mu$. Our method
is based on a companion linearization derived from an accurate Chebyshev
interpolation of $A(\mu)$ on the interval $[-a,a]$, $a \in \mathbb{R}$. The
solution to the linearization is approximated in a preconditioned BiCG setting
for shifted systems, where the Krylov basis matrix is formed once. This process
leads to a short-term recurrence method, where one execution of the algorithm
produces the approximation to $x(\mu)$ for many different values of the
parameter $\mu \in [-a,a]$ simultaneously. In particular, this work proposes
one algorithm which applies a shift-and-invert preconditioner exactly as well
as an algorithm which applies the preconditioner inexactly. The competitiveness
of the algorithms are illustrated with large-scale problems arising from a
finite element discretization of a Helmholtz equation with parameterized
material coefficient. The software used in the simulations is publicly
available online, and thus all our experiments are reproducible.
