Micro- and nanostructuring of optical media is the cornerstone for a wide range of applications in optics. Artificial structures significantly enrich the interaction of light with dielectrics and metals and are harnessed to reveal exciting optical effects. Recently hypercrystals were proposed as intriguing class of periodic structures bringing together hyperbolic dispersion law of hyperbolic metamaterials and the Bragg diffraction. In this work we proposed the finite-difference frequency-domain method for modeling of light diffraction in these highly anisotropic structures. The method is suitable for one- and two-dimensional gratings as well as multilayer gratings made of anisotropic materials. Developed approach gives advantages over the finite-difference time-domain method (FDTD) and rigorous-coupled mode analysis (RCWA) widely used for calculation of light diffraction in gratings. Proposed method eliminates evaluation of time dynamics and discretization of whole volume of a structure as the FDTD does. Instead, Yee’s mesh of the cross section of the unit cell of grating is involved. As well information about the Bloch modes of the field inside the structure inaccessible directly for the FDTD is obtained. We used of values of electromagnetic field in mesh nodes instead Fourier series applied in the RCWA to overcome difficulties of the RCWA for high contrast and highly anisotropic structures as hypercrystals. Using the proposed numerical approach, we demonstrate that interface of one-dimensional hypercrystals formed by ridges of hyperbolic metamaterial can support excitation of the Dyakonov plasmons, which merge the Dyakonov waves supported by the interface of birefringent crystal and surface plasmon polaritons in metals.