In this work we consider the problem of a matrix equation whose matrix and right hand side depend on a set of parameters. We propose a interpolation method that approximates the value of the output at given parameter values by a linear combination of solutions already computed at arbitrarily specified points in the parameter space. To compute the interpolant, we solve an optimization problem that yields the coefficients of the linear combination. The framework is very general and begets a straightforward adaptive strategy for choosing new solutions to add to the approximation based on the largest minimum value of the objective function over the parameter space. Our interpolation based model reduction approach also easily generalizes to nonlinear systems of equations. The efficiency of our approach will be demonstrated with numerical examples on both linear and nonlinear problems.