A Wavelet-Galerkin method is proposed to solve the singular perturbation problem with boundary layers numerically. Because there are boundary layers in the solution of the singular perturbation problem, the approximation spaces with different scale wavelets and boundary bases are chosen. In addition, the computation of the inner integrals is transformed to an eigenvalue problem. Therefore, a high accuracy method with reasonable computation is obtained. On the other hand, there is an explicit diagonal preconditioning which makes the condition number of the stiff matrix become bounded by a constant. The error estimate of the Wavelet-Galerkin method and the analysis of the computation complexity are given. The numerical examples show that the method is feasible and effective for solving the singular perturbation problem with boundary layers numerically.