We present an approach of combining the plane-wave method and the real-space treatment to characterise the periodic variation and the decay of wave functions in the material plane and from the material surfaces. The proposed approach is natural for 2D material systems, and thus may avoid some intrinsic limits involving the artificial replication of material layer in traditional supercell methods. In particular, we show that the proposed method is easy in the implementation and, especially, efficient in the computation since low-cost computational algorithms, such as iterative and recursive techniques, can be invoked to treat matrices with the block tridiagonal structure. Using the proposed approach we show first-principles evidences to supplement the current knowledge of some fundamental issues in bilayer graphene systems, including the coupling between the two graphene layers, the preservation of the σ-band of monolayer graphene in the electronic structure of the bilayer system, and the difference between the low-energy band structure of the bilayer AA- and AB-stacking configurations