Not sure if its too late for this answer, but just in case.
The density of states is defined as $D(E)=\int_{1st BZ}\delta(E-E(\mathbf{k}) d\mathbf{k}$.
For analytic results the following paper was the first one to calculate the DOS in a honeycomb structure with nearest-neighbor hopping (whose dispersion relation is the one that you mention): http://journals.aps.org/pr/abstract/10.1103/PhysRev.89.662 .
For numerical results just apply the definition of $D(E)$. Explicitly you should create two arrays containing values of $k_x$ and $k_y$ (say 100 for each going from their minimum values to their maximum values). Calculate its associated energies. Cut the 1st BZ into same-sized $\mathbf{k}$ regions. Finally calculate the number of states in each area (sum 1 for every $E(\mathbf{k})$ lying inside a given area bounded by the values of $E(\mathbf{k})$ evaluated at the boundaries of each of the $\mathbf{k}$ regions in the 1st BZ). Obviously, the dimensions of your $\mathbf{k}$ area should be greater than the difference between consecutive $k_x$ and $k_y$. Finally normalize your result.