We consider approximating solutions to parameterized linear systems of the form A(μ1,μ2)x(μ1,μ2)=b, where (μ1,μ2)∈R2. Here the matrix A(μ1,μ2)∈Rn×n is nonsingular, large, and sparse and depends nonlinearly on the parameters μ1 and μ2. Specifically, the system arises from a discretization of a partial differential equation and x(μ1,μ2)∈Rn, b∈Rn. The treatment of linear systems with nonlinear dependence on a single parameter has been well-studied, and robust methods combining companion linearization, Krylov subspace iterative methods, and Chebyshev interpolation have enabled fast solution for multiple values of a single parameter at the cost of a single iteration. For systems depending nonlinearly on multiple parameters, the computation becomes much more challenging. This work overcomes those additional challenges by combining robust strategies for treating a single parameter with tensor-base approaches; i.e., combining companion linearization, the Krylov subspace method preconditioned bi-conjugate gradient (BiCG), and a decomposition of a tensor matrix of precomputed solutions, called snapshots. This produces a reduced order model of x(μ1,μ2), and this model can be evaluated inexpensively for many values of the parameters. An interpolation of the model is used to produce approximations on the entire parameter space. In addition this method can be used to solve a parameter estimation problem. This approach allows us to achieve similar computational savings as for the one-parameter case; we can solve for many parameter pairs at the cost of many fewer applications of an efficient iterative method. The technique is presented for dependence on two parameters, but the strategy can be extended to more parameters using the same theoretical approach. Numerical examples of a parameterized Helmholtz equation show the competitiveness of our approach. The simulations are reproducible, and the software is available online.