The initial system is overdetermined: you could simply write
$$\vec D = \alpha \vec E + \beta \vec H \\ \vec B = \gamma \vec E + \delta \vec H,$$
and it describes exactly the same class of fields, but with only four constants instead of the five in your description. As such, you're perfectly free to define $\epsilon=\alpha$, $\mu=\delta$, $c=1/\sqrt{\mu\epsilon}=1/\sqrt{\alpha\delta}$, $\lambda=\beta c=\beta/\sqrt{\alpha\delta}$ and $v=\gamma c=\gamma/\sqrt{\alpha\delta}$, or in other words, you can impose the condition $\mu\epsilon=1/c^2$ for free.
However, it's much less clear to me that you can impose a specific relation between the two mixed coefficients $\lambda$ and $v$, and in general if the material is not reciprocal (cf. this Wikipedia page) then you won't be able to relate the two. Either way, your set-piece explicitly disavows the ability to impose constraints on the five constants you're given, so setting $\lambda=v$ is at odds with the problem as set.
Thus, the best you can do is to take the wave equations as you've derived them,
\begin{align}
-\nabla^2 \hat E &= [\mu \epsilon - (v/c)^2] \omega^2 \hat E + \mu / c [\lambda - v] \omega^2 \hat H \\
\nabla^2 \hat H &= \epsilon / c [\lambda - v] \omega^2 \hat E + [(\lambda/c)^2 - \mu \epsilon] \omega^2 \hat H,
\end{align}
and then assume that both fields are in Helmholtz eigenstates of the laplacian, giving you
\begin{align}
k^2 \hat E &= [\mu \epsilon - (v/c)^2] \omega^2 \hat E + \mu / c [\lambda - v] \omega^2 \hat H \\
-k^2 \hat H &= \epsilon / c [\lambda - v] \omega^2 \hat E + [(\lambda/c)^2 - \mu \epsilon] \omega^2 \hat H
\end{align}
or in other words
$$
\begin{pmatrix}
k^2-\left(\mu\epsilon-\frac{v^2}{c^2}\right)\omega^2 &
\frac{\mu}{c}(\lambda -v)\omega^2 \\
\frac{\epsilon}{c}(\lambda -v)\omega^2 &
k^2-\left(\mu\epsilon-\frac{\lambda^2}{c^2}\right)\omega^2
\end{pmatrix}
\begin{pmatrix}
\vec E\\\vec H
\end{pmatrix}
=0,
$$
which then requires that the determinant vanish, i.e.
$$
\left(k^2-\left(\mu\epsilon-\frac{v^2}{c^2}\right)\omega^2 \right)
\left(k^2-\left(\mu\epsilon-\frac{\lambda^2}{c^2}\right)\omega^2 \right)
-
\frac{\mu}{c}(\lambda -v)\omega^2
\times\frac{\epsilon}{c}(\lambda -v)\omega^2
=0
$$
which fails to simplify in ways that make suspect you've made a sign error somewhere, so I'll leave it to you to double-check your workings and pull out the dispersion relation from there.
