An efficient numerical method is developed for calculating surface modes of a semi-infinite two dimensional photonic crystals by Dirichlet-to-Neumann (DtN) maps.Using the DtN maps of three different unit cells to construct the DtN map of the supercell,the surface modes are solved from a linear eigenvalue problem on two boundaries of the supercell,leading to relatively small matrices.Finally we use two numerical examples to verify the effectiveness of the method.