We propose a fast and accurate numerical method to price European swaptions in a multi-factor Gaussian term structure model that can speed up the calibration to the volatility surface. Pricing an interest rate option in such a multi-factor term structure model involves evaluating a multi-dimensional integration of the underlying claim's payo on a domain where the payo is positive. In our method, we approximate the exercise boundary of the state variables as a hyperplane tangent to the maximum probability point on the boundary and simplify the multi-dimensional integration into an analytic form. The maximum probability point can be found using the gradient descent method. We show that the quality of our method is superior to the result of previous studies by directly comparing them to the exact price obtained from the numerical integration.