I don't understand your question. Once you solve the eigenvalue problem ##D - \omega^2 I = 0## (see the previous answer), you find four frequencies as a function of your wave-vector ##\mathbf k##:
$$\omega_i = \omega_i(\mathbf k) = \omega_i(k_x, k_y)$$
since this a two-dimensional problem, i.e. ##\mathbf k = (k_x, k_y)##. Now, if you are interested in the frequencies at a particular point of the Brillouin zone, you just plug in the coordinates of that point: for example, if you want to calculate the frequencies at M, which has coordinates ##(2 \pi/\sqrt{3}a, 0)##, you just evaluate ##\omega_i(2 \pi/\sqrt{3}a, 0)##.
If you want to plot the dispersion along the direction ##\Gamma \rightarrow M##, since both ##\Gamma## and ##M## have ##k_y = 0##, you simply need to calculate ##\omega_i(k_x, 0)## and, as in the paper, you'll find 4 mathematical expression as functions of ##k_x## (see eq. 26-29).