In this paper we present a numerically stable method for the model order reduction of finite element (FE) approximations to passive microwave structures parameterized by polynomials in several variables. The proposed method is a projection-based approach using Krylov subspaces and extends the works of Gunupudi et al. (P. Gunupudi, R. Khazaka and M. Nakhla, Analysis of transmission line circuits using multidimensional model reduction techniques, IEEE Trans. Adv. Packaging 25 (2002), pp. 174-180) and Slone et al. (R.D. Slone, R. Lee and J.-F. Lee, Broadband model order reduction of polynomial matrix equations using single-point well-conditioned asymptotic waveform evaluation: derivations and theory, Int. J. Numer. Meth. Eng. 58 (2003), pp. 2325-2342). First, we present the multivariate Krylov space of higher order associated with a parameter-dependent righthand-side vector and derive a general recursion for generating its basis. Next, we propose an advanced algorithm to compute such a basis in a numerically stable way. Finally, we apply the Krylov basis to construct a reduced order model of the moment-matching type. The resulting single-point method requires one matrix factorization only. Numerical examples demonstrate the efficiency and reliability of our approach.