How to build numerically a dispersion relation of a finite phononic system?

I have doubts about how to build the plot of dispersion relation for a general finite phononic system, in particular spring-mass systems. Eqs of any of these systems can be written in matrix form as $$D\vec{u}=\omega^2 \vec{u} \tag{1}$$ For the case of "infinite" system, usually one set the periodic boundary condition (PBC) and under the Bloch expansion $$u_i=u_0 \mathrm{e}^{ikn}$$ in (1), one can solve analytically $$\omega(k)$$ for a single primitive cell. As example, in the 1d diatomic chain one can obtain $$\omega(k)$$ by this way

In this case, considering $$N$$ masses and the PBC, this quantizes the values of $$k$$, yielding $$N$$ points in $$\omega(k)$$, but taking $$N\to\infty$$, one can plot a continuous line in $$\omega(k)$$ (similar for monoatomic chain).

However in a finite system, under specific boundary conditions in $$D$$, if I impose the Bloch expansion in (1), the lack of translational invariance of $$D$$ doesn't allow us to reduce (1) to the dynamics of a single cell. I'm not interested in analytically solving $$\omega(k)$$ in these system, but I would like to be able to plot the N points of $$\omega(k)$$. I think that I'm wrong imposing $$u_0 \mathrm{e}^{ikn}$$ here, but then I don't know how to make appear the wave vector $$k$$ in (1) to obtain the $$\omega(k)$$ relations. Is there a generic numerical way to obtain the points of $$\omega(k)$$ in these systems?

(e.g. the diatomic chain with one end free and the other fixed)

You can still use periodic boundary conditions for a finite chain if the ends are fixed (as I will demonstrate below), and so the dispersion relation is unchanged in the sense that you just replace the continuous variable $$k$$ with its discrete version.
To treat a finite chain first impose periodicity as you did before, this means that the two ends of the chain $$x_0$$ and $$x_N$$ will be regarded as the same point but this is okay since we fix both ends by demanding that $$x_0(t) = x_N(t) = 0$$ If we look at the general solution $$x_n(t) = \sum_{k=-\pi/a}^{\pi/a}e^{ikna}(A_ke^{i\Omega_k t}+B_ke^{-i\Omega_k t})$$where $$\Omega_k$$ is the mode frequency, the ends are fixed if we set $$A_k=-A_{-k}$$, $$~B_k=-B_{-k}$$, and $$A_0=B_0=0$$.