Current methods of parametric model order reduction based on the interpolation of system matrices are extended in this paper to an efficient method to tackle multidimensional parameter spaces. In order to overcome the curse of dimensionality the technique of sparse grids with multiple outputs is applied. The procedure is divided into an offline and online phase. In the offline phase the weighting matrices to the corresponding basis functions are computed with respect to generalized coordinates and appropriate matrix manifolds. The number of levels for the sparse grid can be set manually or determined by a tolerance criterion. The calculation of the interpolated system takes place during the online phase by evaluating the interpolant. The performance of the proposed method is demonstrated by two numerical examples.