I agree with you completely, in fact I came to a similar conclusion. I think that there is an implicit (or perhaps explicit somewhere I overlooked) assumption when solving linear homogeneous systems namely that the unknowns are somehow "independent" in the sense that picking the value of one does not automatically pick the value of the other. In the case of general quadratic systems involving mixed terms (not the simple case of a diagnal matrix I inquired about initially) it is possible in principle to replace each (degree 2) mixed term with a new variable, resulting in total in a new system of n (n - 1) / 2 variables. However whereas in my simple original case the domain of possible solutions is R^n, in the new system I just constructed the sets of values the variables can take on independent of the coefficient matrix [s_ij] and the right hand side constraint of the system (call it B) is just some (possibly disconnected?) subset of R^n (i.e. because of their interdependence). So the task of find the unique solution w.r.t the coefficient matrix [s_ij] and B involves searching through a subset of R^n with a rather complex (and exponential) structure which definitely cannot be done by a Gaussian Elimination.
But that is a very loose reasoning I am giving. Not even sure my reasoning was even correct. In any case I will follow up your suggestion and play a bit