In order to find the solutions of the system of conics you mention, it suffices to give a procedure to find the projection of the simultaneous vanishing set of the two conics onto some axis, for instance, the $y$-axis: the coordinates of the projection onto the $y$-axis are the $y$-coordinates of the intersection points. Knowing them allows you to substitute back into the equations the values of $y$ and solve a system of equations in a single variable $x$, thus making the problem simpler. In terms of ideals, suppose that we can find a non-zero element $r$ in the ideal generated by the two conics, that depends only on a single variable, say $y$. This means that every solution of the system has $y$-coordinate satisfying the polynomial $r$. Thus we have severely limited the choices for the $y$-coordinates of the intersection (namely, they all have to satisfy the polynomial $r$) and, provided we can actually solve the polynomial $r$ in one variable, we can then substitute the various values of $y$ we found back into the initial equations and solve for $x$, again using our algorithm for solving polynomials in a single variable. So far, so good, I hope! The question is how to produce the element $r$. This is where Gröbner bases come in.
In order to do the computation you ask, one would have to choose an appropriate monomial order. I will gloss over this, and simply "do the obvious". It is an exercise for you to figure out which order to use so that the computation I am going to make is actually a Gröbner bases computation. Scattered throughout the computation there will be also a couple of special cases that I will not deal with: again, you can treat those as exercises for you!
First, I am going to choose a generic basis, so that the first conic has the form
$$
x^2 + \alpha x + \beta ,
$$
where $\alpha$ and $\beta$ are polynomials in $y$ only (I am trying to simplify the notation as much as possible; the only "assumption that I have made is that the coefficient of $x^2$ is non-zero, which can be arranged unless the "conic" is defined by a polynomial of degree at most one).
In this basis, the second equation can be assumed to have no $x^2$ term (indeed, eliminating the $x^2$ term using the first equation would be the first step in any reasonable Gröbner basis computation in which the term $x^2$ is the highest term in sight). Thus I am going to write the second conic as
$$
\gamma x + \delta
$$
where, as before, $\gamma$ and $\delta$ are polynomials in $y$ only. Again, just to fix on a definite case, I am going to assume that the polynomial $\gamma$ is non-zero. (If $\gamma$ were zero, then we would have found a polynomial in the ideal generated by the two conics which only depends on $y$: this was our goal at the start anyway! Obviously, you should worry about the case $\gamma=\delta=0$, but I won't.)
To compute a Gröbner basis, you would now compute S-polynomials: let's do it here as well. I am going to try to eliminate the $x^2$ term from the first equation, by using the second equation. This is easy: multiply the first equation by $\gamma$ and the second one by $x$ and take the difference: we are left with the equation
$$
(\alpha \gamma - \delta) x + \beta \gamma .
$$
Now we are going to eliminate the $x$ term from this last equation using again the second equation: multiply the second equation by $(\alpha \gamma - \delta)$, the last equation by $\gamma$ and subtract to obtain
$$
\delta^2 - \alpha \gamma \delta + \beta \gamma^2 .
$$
We found an expression independent of $x$!! We are done... provided this expression is not identically zero. You can figure out what this would mean and what happens in this case. Note also that the final expression is what we would have obtained if, at the very beginning, we had "solved" $x=-\frac{\delta}{\gamma}$ using the second equation, substituted in the first equation and cleared the denominators.
I hope that this "hybrid" computation explains what is going on: Gröbner bases and Buchberger's algorithm are a systematic way of "solving" systems of equations. You do not have to do any thinking, once you set up the problem. But you need to set up the problem so that it computes what you want. In this case, you could have used several shortcuts to get to the answer, without following all the steps. In more complicated situations, Buchberger's algorithm might be the best way of keeping track of all the steps to be taken.
Let me also comment that, except in the case in which you have a computer doing the computations for you, it is highly unlikely that Gröbner bases will help you with a specific question, unless you could have also found simple tricks to solve it right away.
If $\alpha$ is a root where the graph of $f$ is tangent to the $x$ axis, then $\alpha$ is a root of multiplicity $\ge 2$ (not necessarily even).
Case 1: $\alpha$ has multiplicity $3$. Then $\alpha$ is a root of the polynomial $f''(x)$ which has degree $1$ and rational coefficients, and therefore $\alpha$ is rational.
Case 2: $\alpha$ has multiplicity $2$. Then the gcd of $f(x)$ and $f'(x)$ is $x-\alpha$. But the gcd of two polynomials with rational coefficients has rational coefficients, so $\alpha$ must be rational.
Best Answer
The problem is actually fairly simple if the operations are performed in the right order. When $gg_i$ and the $g_jf_x$ polynomials are taken modulo $f$ first, the remainders' monomials will have degrees of the same order (
rem(f,g)in SymPy), so that a linear system can be set up to find the $a_{ij}$. To illustrate, for the example's polynomials, the entries of $A$'s first row are the solutions to the linear system starting with $$\begin{bmatrix} 0&12&-4\\ 106&-108&196\\ 96&-128&200\\ \vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}a_{11}\\a_{12}\\a_{13}\end{bmatrix}=\begin{bmatrix}0\\4\\6\\\vdots\end{bmatrix}$$ where the columns from left to right are the reduced $g_jf_x$ and $gg_1$ polynomials respectively, and the displayed rows correspond to the $x^4y^3,x^4y^2,x^4y$ coefficients. Once that hurdle was cleared, I managed to complete the implementation and successfully find splitting fields for small bivariate polynomials.Immediately afterwards, however, I saw an improvement from Jürgen Gerhard mentioned in the same paper that saves the hassle of finding a basis for the initial linear system $G$'s nullspace, guessing $g$ and constructing $A$ – the hassle that led me to ask this question in the first place. It involves taking any non-trivial $g$ in $G$ and computing the resultant $\operatorname{Res}_x(f,g-zf_x)$, from which the number of absolutely irreducible factors and the splitting field can be derived. The sheer size of the $G$ matrices I encountered with slightly larger polynomials also compelled me to write the implementation in PARI/GP instead, which managed to find a quintic splitting field for a degree-$25$ test polynomial I concocted. The PARI/GP implementation is available as
gao.gphere.