First, note that we are dealing with polynomial functions here so we can freely impose nonvanishing of certain expressions (as long as there is some $\vec{a},\vec{b},\vec{c}$ that achieves this) and then use a continuity argument to pass to the limiting case. You have already use this when you say you can write the triple product as a linear combination of $\vec{b}$ and $\vec{c}$ (you have thrown away cases where $\vec{b}\parallel\vec{c}$ only to regain them at the limit).
To start with, assume $\vec{a}\cdot\vec{b},\vec{a}\cdot\vec{c}$ are not simultaneously zero. Then the solutions $(x,y)$ to $0=(\vec{a}\cdot\vec{b})x+(\vec{a}\cdot\vec{c})y$ lies on the line perpendicular to $(\vec{a}\cdot\vec{b},\vec{a}\cdot\vec{c})$ and hence you can introduce a constant of proportionality $\lambda$ (a priori could depend on $\vec{a},\vec{b},\vec{c}$) so that
$$
(x,y)=\lambda(\vec{a}\cdot\vec{c},-\vec{a}\cdot\vec{b})
$$
and what's more, this constant $\lambda$ can be chosen to be independent of $\vec{a},\vec{b},\vec{c}$. Your reference just assert this and went on to calculate $\lambda$ using specific case. But here is a simpler way to justify part of it that would suffice for the derivation.
We can further restrict our attention to the case $\vec{a}\times(\vec{b}\times\vec{c})\neq\vec{0}$. Then it doesn't hurt to replace $\vec{a}$ by its orthogonal projection along $\vec{b}\times\vec{c}$ to the plane containing $\vec{b}$ and $\vec{c}$ (because it wouldn't change the vector triple product nor $\lambda$).
Hence we just need to show that $\lambda$ is the same for $\vec{b}\times(\vec{b}\times\vec{c})$ and $\vec{c}\times(\vec{b}\times\vec{c})$. Taking dot product with $\vec{c}$ in the $\vec{b}\times(\vec{b}\times\vec{c})$ gives
$$
[\vec{c},\vec{b},\vec{b}\times\vec{c}]=
\lambda(\vec{c}\cdot\vec{b})(\vec{b}\cdot\vec{c})-\lambda(\vec{c}\cdot\vec{c})(\vec{b}\cdot\vec{b})
$$
but the LHS is, by cyclic permutation of the scalar triple product
$$
\begin{align*}
\vec{c}\cdot(\vec{b}\times(\vec{b}\times\vec{c}))
&=(\vec{b}\times\vec{c})\cdot(\vec{c}\times\vec{b})\\
&=-|\vec{b}\times\vec{c}|^2\\
&=-b^2c^2\sin^2\theta\\
&=(bc\cos\theta)^2-b^2c^2\\
&=(\vec{b}\cdot\vec{c})^2-(\vec{b}\cdot\vec{b})(\vec{c}\cdot\vec{c})
\end{align*}
$$
So $\lambda=1$ in this case.
Similarly, taking scalar product with $\vec{b}$ for $\vec{c}\times(\vec{b}\times\vec{c})$, we get $\lambda=1$ too in this case.
So we can take $\lambda=1$ for all $\vec{a},\vec{b},\vec{c}$.