Many engineered photonic devices can be decomposed into parts where the material properties are independent of one or more spatial variables. Numerical mode-matching methods are widely used to simulate such photonic devices due to the efficiency gained by treating the separated variables analytically. Existing mode-matching methods based on piecewise polynomials are more accurate than those based on the global Fourier basis or low-order finite difference, finite-element schemes, but they may exhibit numerical instability when a large number of eigenmodes are used. To overcome this difficulty, we introduce the spectral Galerkin mode matching method (SGMM) based on a global piecewise-polynomial basis and a Galerkin method to solve the eigenmodes. It is shown that the numerical eigenmodes of SGMM preserve the pseudo-orthogonality of the analytical eigenmodes. This property leads to linear systems that are typically well-conditioned. Numerical examples indicate that SGMM is more stable than other mode matching methods, and gives reliable results even when a large number of eigenmodes are used.