In this paper (open access), the authors used Fourier series with most general wave equation to find the dispersion relation. I am presenting some main equations as snippets to depict their solution. We want to solve for $u(x,t)$ in the following wave equation (Eq. (1) in the paper):
\begin{equation}\tag{1} \frac{\partial}{\partial x}\left[ E(x,t)\frac{\partial}{\partial x} u(x,t) \right]-\frac{\partial}{\partial t}\left[ \rho(x,t)\frac{\partial}{\partial t} u(x,t) \right]=0. \end{equation}
So we expand or represent the quantities using Fourier series:
$$\tag{2} E(x,t)=\sum_{p=-\infty}^\infty E_{p}e^{ip(\omega_mt-k_mx)}, $$ and $$\tag{3} \rho(x,t)=\sum_{p=-\infty}^\infty \rho_{p}e^{ip(\omega_mt-k_mx)}, $$ both having periodicty of $\omega_m$ and $k_m$ for time and space. Now we assume a general solution (Floquet form) for $u(x,t)$ as:
$$\tag{4} u(x,t)=e^{i(\omega t-kx)} \sum_{n=-\infty}^\infty u_{n}e^{in(\omega_mt-k_mx)}, $$
and plug everything in the wave-equation. The authors say that using orthogonality we directly have the equation (9 in the paper):
$$\tag{5} \sum_{n=-\infty}^\infty [(k+nk_m)(k+pk_m)]E_{p-n}u_n = \sum_{n=-\infty}^\infty [(\omega+n\omega_m)(\omega+p\omega_m)]\rho_{p-n}u_n.$$
Then they say that it can be "recast" in the form (Eq. 11 in paper):
$$\tag{6} (\omega^2 \hat{L}_2 + \omega \hat{L}_1 + \hat{L}_0)\vec{u}=0, $$ which is a matrix equation in $\vec{u}$ and connects $u_n$, $u_{n−1}$, and $u_{n+1}$.
How can Eq. 5 be cast in Eq. 6 given that Eq. 5 does not have a relation between $u_n$, $u_{n−1}$, and $u_{n+1}$? It only has two independent sums on both sides. They do not have different phases on both side to shift a phase term and take common summation to relate these.
