TL;DR
With Matlab or Mathematika I cannot help, but
in Python there is an implementation available: sc_dos
Here $D$ is the half-bandwidth $D = 6t$.
import numpy as np
import gftool as gt
eps = np.linspace(-1.2, 1.2, num=6001)
dos = gt.sc_dos(eps, half_bandwidth=1)
The evaluation of the DOS takes for me ~100 ms.
You already gave the correct expressions. We have the Green's function
$$
G(z) = \frac{1}{N} \sum_{\boldsymbol{k}} \frac{1}{z-\epsilon_{\boldsymbol{k}}}
$$
and the (normalized) density of states (DOS)
$$
D(\epsilon)
= \frac{1}{N} \sum_{\boldsymbol{k}} \delta(\epsilon - \epsilon_{\boldsymbol{k}})
= -\frac{1}{\pi} \Im G(\epsilon+i0^+),
$$
where $\epsilon$ is a real energy variable.
The second equality is Sokhotski–Plemelj.
A naive sum over the points is extremely demanding, as a tremendous number of $\boldsymbol{k}$ points is necessary in 3 Dimensions.
To smoothen the function, we can evaluate the Green's function on a contour parallel to the real axis shifted by a finite $\eta>0$ into to upper complex half-plane:
$$
D(\epsilon) \approx -\frac{1}{\pi} \Im G(\epsilon+i\eta).
$$
The bigger we choose $\eta$ the smoother the function becomes, but on the other hand we loose features.
As we are only interested in the thermodynamic limit $N\rightarrow \infty$,
a smarter approach than just sampling $\boldsymbol{k}$, is to replace the sum by the integral. For integrals, we have more or less efficient algorithms.
So let's calculate
$$
G(z) = \int \frac{d^3 k}{{2\pi}^3} \frac{1}{z-\epsilon_{\boldsymbol{k}}}
$$
instead. The $\epsilon_\boldsymbol{k}$, is symmetric for all $k_{x_i}$: $\epsilon(k_{x_i}) = \epsilon(-k_{x_i})$, thus it is enough to integrate over a eighth of the Brillouin zone.
And finally we can use analytic results for the integrals. We note that we can express the 3D Green's function in terms of known results of the 1D and 2D Green's function as we have
$$
\epsilon_{\boldsymbol{k}} = \epsilon_{k_x, k_y, k_z}
= \epsilon^{2D}_{k_x, k_y} - 2t \cos(k_z)
= \epsilon^{1D}_{k_x} - 2t[\cos(k_y) + \cos(k_z)]
$$
and therefore
$$
G(z)
= \int \frac{dk_z}{2\pi} G^{2D} (z - 2t\cos(k_z))
= \int \frac{dk_y}{2\pi} \int \frac{dk_z}{2\pi} G^{1D}(z - 2t[\cos(k_y) + \cos(k_z)]).
$$
The one-dimensional Green's function $G^{1D}(z)$ can be easily evaluated, the two-dimensional Green's function $G^{2D}(z)$ can be expressed in terms of the complete elliptic integral of first kind (which can be found in standard text-books). Using $G^{2D}(z)$ is basically the result given by bRost03.
A very smart guy named Joyce even found an expression for $G^{3D}(z)$ in 1973.
The equations are a bit lengthy and complicated, so I will avoid copying them here. But we implemented them in a Python module gftool>=0.8.0, see sc_dos. You will also find the relevant references there.
