Abstract. We propose an asymptotic numerical method called the Gaussian beam approach for the boundary value problem of high frequency Helmholtz equation. The basic idea is to approximate the traveling waves with a summation of Gaussian beams by the least squares algorithm. Gaussian beams are asymptotic solutions of linear wave equations in the high frequency regime. We deduce the ODE systems satisfied by the Gaussian beams up to third order. The key ingredient of the proposed method is the construction of a finite-dimensional beam space which has a good approximating property. If the exact solutions of boundary value problems contain some strongly evanescent wave modes, the Gaussian beam approach might fail. To remedy this problem, we resort to the domain decomposition technique to separate the domain of definition into a boundary layer region and its complementary interior region. The former is handled by a domain-based discretization method, and the latter by the Gaussian beam approach. Schwarz iterations should then be performed based on suitable transmission boundary conditions at the interface of two regions. Numerical tests demonstrate that the proposed method is very promising.