using this subroutine over other methods, for instance the Newton’s method. This subroutine does not require initial guesses for the roots but only the contour which encloses all the roots of the polynomial. Besides, it also avoids calculating the derivatives which may result into numerical overflow.

In alternative to [*ezero*]{} we have used a different approach for calculating the roots. We survey the complex $\tilde{k}$-plane for the zeros of the $\text{Det}[M(\tilde{k})M(\tilde{k})^{\dagger}]$ where $M_{rs}=-2\cos(\tilde{k})\delta_{rs}-\tilde{\mathcal{H}}_{rs}$ for $r,s=1,...,N$ [@Neuberger]. (In our system $-\pi<\Re\{\tilde{k}\}<\pi$ and $\Im\{\tilde{k}\}<0$.) These zeros give the eigenvalues of $\tilde{\mathcal{H}}$. However, in the latter approach it is advisable to disintegrate the complex plane into small cells at first and then at every iteration into smaller one - only for $N$ cells which contain minima of the lowest eigenvalue and throwing the rest out. In this way one makes the algorithm faster and obtain the zeros in a reasonable precision. For a tridiagonal matrix this algorithm consumes a time which roughly grows with $N^3$. However, while comparing the two methods on a simple machine we find that the method used in [*ezero*]{} is much faster than the method described here. We refer to [@ezero] for further details of this subroutine.

In Fig. 