Summary: We consider linear parameterized systems \(A(\mu)x(\mu)=b\) for many different \(\mu\), where \(A\) is large and sparse and depends nonlinearly on \(\mu\). Solving such systems individually for each \(\mu\) would require great computational effort. In this work we propose to compute a partial parameterization \(\tilde{x}\approx x(\mu)\), where \(\tilde{x}(\mu)\) is cheap to evaluate for many \(\mu\). Our methods are based on the observation that a companion linearization can be formed where the dependence on \(\mu\) is only linear. In particular, methods are presented that combine the well-established Krylov subspace method for linear systems, GMRES, with algorithms for nonlinear eigenvalue problems (NEPs) to generate a basis for the Krylov subspace. Within this new approach, the basis matrix is constructed in three different ways, using a tensor structure and exploiting that certain problems have low-rank properties. The methods are analyzed analogously to the standard convergence theory for the method GMRES for linear systems. More specifically, the error is estimated based on the magnitude of the parameter \(\mu\) and the spectrum of the linear companion matrix, which corresponds to the reciprocal solutions to the corresponding NEP. Numerical experiments illustrate the competitiveness of the methods for large-scale problems. The simulations are reproducible and publicly available online.