Determining orthonormal eigenvectors of the DFT matrix, which is closer to the samples of Hermite-Gaussian functions, is crucial in the definition of the discrete fractional Fourier transform. In this work, we disclose eigenvectors of the DFT matrix inspired by the ideas behind bilinear transform. The bilinear transform maps the analog space to the discrete sample space. As j omega in the analog s-domain is mapped to the unit circle one-to-one without aliasing in the discrete z-domain, it is appropriate to use it in the discretization of the eigenfunctions of the Fourier transform. We obtain Hermite-Gaussian-like eigenvectors of the DFT matrix. For this purpose we propose three different methods and analyze their stability conditions. These methods include better conditioned commuting matrices and higher order methods. We confirm the results with extensive simulations.