It has been difficult to solve eigenmodes of plasmonic crystals in two or three dimensions either analytically or numerically. In this study, we present an interfacial opera- tor approach for solving guided wave modes of plasmonic crystals. They are formulated as an eigenvalue problem of the wavenumber along the axis of the crystal. In this formulation, the permittivity and permeability of the metallic component can be arbitrary functions of frequency. Moreover, a coupling interface method is introduced to facilitate accurate treatment of the in- terface conditions with an arbitrary shape between the metal and host materials. Numerical results are illustrated for different shapes of plasmonic crystals, layered, cylindrical and split-ring structures. The physical significance is discussed. Finally, it is demonstrated that the present method can resolve fine eigenmodes of the split-ring structure.