Global linear instability analysis is a powerful tool for the complex flow diagnosis. However, the methods used in the past would generally suffer from some disadvantages, either the excessive computational resources for the low-order methods or the tedious mathematical derivations for the high-order methods. The present work proposed a CFD-aided Galerkin methodology which combines the merits from both the low-order and high-order methods, where the expansion on proper basis functions is preserved to ensure a small matrix size, while the differentials, incompressibility constraints and boundary conditions are realized by applying the low-order linearized Navier-Stokes equation solvers on the basis functions on a fine grid. Several test cases have shown that the new method can get satisfactory results for one-dimensional, two-dimensional and three-dimensional flow problems and also for the problems with complex geometries and boundary conditions.