We propose a fast and accurate numerical method for pricing European swaptions in multifactor Gaussian term structure models. Our method can be used to accelerate the calibration of such models to the volatility surface. The pricing of an interest rate option in such a model involves evaluating a multidimensional integral of the payoff of the claim on a domain where the payoff is positive. In our method, we approximate the exercise boundary of the state space by a hyperplane tangent to the maximum probability point on the boundary and simplify the multidimensional integration into an analytical form. The maximum probability point can be determined using the gradient descent method. We demonstrate that our method is superior to previous methods by comparing the results to the price obtained by numerical integration.