In this study, we introduce an efficient method for determining isolated singular points of two-dimensional semi-infinite and bi-infinite photonic crystals, equipped with perfect electric conductor and quasi-periodic mixed boundary conditions. This specific problem can be modeled by a Helmholtz equation and is recast as a generalized eigenvalue problem involving an infinite-dimensional block quasi-Toeplitz matrix. Through an intelligent implementation of cyclic structure-preserving matrix transformations, the contour integral method is elegantly employed to calculate the isolated eigenvalue and to extract a component of the associated eigenvector. Moreover, a propagation formula for electromagnetic fields is derived. This formulation enables rapid computation of field distributions across the expansive semi-infinite and bi-infinite domains, thus highlighting the attributes of edge states. The preliminary MATLAB implementation is available at https://github.com/FAME-GPU/2D_Semi-infinite_PhC.