Dynamical systems are mathematical models characterized by a set of differential or difference equations. Due to the increasing demand for more accuracy, the number of equations involved may reach the order of thousands and even millions. With so many equations, it often becomes computationally cumbersome to work with these large-scale dynamical systems. Model reduction aims to replace the original system with a reduced system of significantly smaller dimension which will still describe the important dynamics of the large-scale model. Interpolation is one method used to obtain the reduced order model. This requires that the reduced order model interpolates the full order model at selected interpolation points. Reduced order models are obtained through the Krylov reduction process, which involves solving a sequence of linear systems. The Iterative Rational Krylov Algorithm (IRKA) iterates this Krylov reduction process to obtain an optimal Η₂ reduced model. Especially in the large-scale setting, these linear systems often require employing inexact solves. The aim of this thesis is to investigate the impact of inexact solves on interpolatory model reduction. We considered preconditioning the linear systems, varying the stopping tolerances, employing GMRES and BiCG as the inexact solvers, and using different initial shift selections. For just one step of Krylov reduction, we verified theoretical properties of the interpolation error. Also, we found a linear improvement in the subspace angles between the inexact and exact subspaces provided that a good shift selection was used. For a poor shift selection, these angles often remained of the same order regardless of how accurately the linear systems were solved. These patterns were reflected in Η₂ and Η∞ errors between the inexact and exact subspaces, since these errors improved linearly with a good shift selection and were typically of the same order with a poor shift. We found that the shift selection also influenced the overall model reduction error between the full model and inexact model as these error norms were often several orders larger when a poor shift selection was used. For a given shift selection, the overall model reduction error typically remained of the same order for tolerances smaller than 1 x 10^-3, which suggests that larger tolerances for the inexact solver may be used without necessarily augmenting the model reduction error. With preconditioned linear systems as well as BiCG, we found smaller errors between the inexact and exact models while the order of the overall model reduction error remained the same. With IRKA, we observed similar patterns as with just one step of Krylov reduction. However, we also found additional benefits associated with using an initial guess in the inexact solve and by varying the tolerance of the inexact solve. / Master of Science