General problem:
If I have some $n \times n$ matrices $\mathsf{M}^\tau$, and column vectors (with $n$ rows) $X^\tau$ is there some mathematical tricks I can do to solve the eigenvalue equation $ \sum_\tau \mathsf{M}^\tau X^\tau = \sum_\tau E^\tau X^\tau$? If there was no summation I would do the usual trick of calculating $\det(\mathsf{M}^\tau - I_n E^\tau) = 0$ but not sure what to do here.
My specific problem:
I have the equation:
$$ \sum_\tau \begin{pmatrix} A+2D & -iC^\tau & -2D & -iC^\tau\\ iC^\tau & B^\tau & -iC^\tau & 0\\ 2D & -iC^\tau & -A-2D & -iC^\tau\\ -iC^\tau & 0 & iC^\tau & -B^\tau \end{pmatrix} \begin{pmatrix} w \\ x^\tau \\ y \\ z^\tau \end{pmatrix} = \sum_\tau E^\tau \begin{pmatrix} w \\ x^\tau \\ y \\ z^\tau \end{pmatrix} $$ and want to find $E^\tau$. $A,B^\tau,C^\tau$ and $D$ are known. The column vectors are unknown.
In these equations $\tau$ is just an index.
Edit: More info:
This is likely superfluous but I'll briefly explain where this equation comes from. Maybe it adds additional information which illuminates what kind of solution can be used.
I have a Hamiltonian in second quantisation representation as
$H = \sum_q \sum_\tau \Bigg\{ A c_q^\dagger c_q + B^\tau \beta_q^{\tau\dagger} \beta_q^\tau + C_q^\tau (c_q + c_q^\dagger )(\beta_{-q}^\tau - \beta_q^{\tau \dagger}) + D(2c_q^\dagger c_q + c_q c_{-q} + c_q^\dagger c_{-q}^\dagger) \Bigg\}$
where the $\beta$'s and $c$'s are creation/annihilation operators. I introduce a new operator $\gamma_q = \sum_\tau w c_q + x^\tau \beta_q^\tau + y c_{-q}^\dagger + z^\tau \beta_{-q}^{\tau \dagger}$ which diagonalises the Hamiltonian if it obeys the Heisenberg equation of motion $[\gamma_q,H]=E^\tau \gamma_q$. Crunching through this equation and representing the result in matrix form gives me the equation above.
I assume that the $X^{\tau}$'s and the $E^{\tau}$'s are both unknown.
One type of solution is found by removing the summation and assuming that $M^{\tau} X^{\tau} = E^{\tau} X^{\tau}$ holds for each individual term. That is, letting $(X^{\tau},E^{\tau})$ be an arbitrary eigenvector-eigenvalue pair of $M^{\tau}$ yields one solution.
Another type of solution is obtained with the Ansatz $X^{\tau} = X$ and $E^{\tau} = E$ (independent of $\tau$). Then it reduces to finding the eigenvectors and eigenvalues of the sum $\sum_{\tau} M^{\tau}$.
I do not know how many other solutions there are. Perhaps it is also possible to pair-wise solutions, e.g. $M^1 X^1 = E^2 X^2$ and $M^2 X^2 = E^1 X^1$, or cyclic solutions $M^k X^k = E^{k+1} X^{k+1}$?
EDIT: Your specific problem looks a little more restricted in that parts of $X^{\tau}$ (namely w and y) are independent of $\tau$. At least the Ansatz $X^{\tau} = X$ should give you one solution here too.