Iterative Greens function, based on cyclic reduction of block tridiagonal matrices, has been the ideal algorithm, through tight-binding models, to compute the surface density-of-states of semi-infinite topological electronic materials. In this paper, we apply this method to photonic and acoustic crystals, using finite-element discretizations and a generalized eigenvalue formulation, to calculate the local density-of-states on a single surface of semi-infinite lattices. The three-dimensional (3D) examples of gapless helicoidal surface states in Weyl and Dirac crystals are shown and the computational cost, convergence and accuracy are analyzed.