I am interested in solving the following system of two quadratic equations in two unknowns, $x$ and $y$,
$$ \begin{cases}A_1x^2 + B_1xy + C_1 y^2 + D_1x + E_1y + F_1 = 0 \\ A_2x^2 + B_2xy + C_2y^2 + D_2x + E_2y + F_2 = 0 \end{cases} $$
Where $A_1, B_1, C_1, D_1, E_1, F_1, A_2, B_2, C_2, D_2, E_2, F_2$ are known given constants.
Is there a systematic algorithm to find all the possible solutions $(x,y)$ ?
You have the following system of quadratic equations:
$ A_1 x^2 + B_1 xy + C_1 y^2 + D_1 x + E_1 y + F_1 = 0$
and
$ A_2 x^2 + B_2 xy + C_2 y^2 + D_2 x + E_2 y + F_2 = 0$
First, transform each of the given quadratic equations into quadratic form notation as follows. Define $ r = \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} $ , then the first equation can be written as the quadratic form $ r^T Q_1 r = 0 $ where
$Q_1 = \begin{bmatrix} A_1 && \dfrac{B_1}{2} && \dfrac{D_1}{2} \\ \dfrac{B_1}{2} && C_1 && \dfrac{E_1}{2} \\ \dfrac{D_1}{2} && \dfrac{E_1}{2} && F_1 \end{bmatrix}$
Similarly, the second quadratic equation can be expressed as the quadratic form $r^T Q_2 r = 0 $, with
$Q_2 = \begin{bmatrix} A_2 && \dfrac{B_2}{2} && \dfrac{D_2}{2} \\ \dfrac{B_2}{2} && C_2 && \dfrac{E_2}{2} \\ \dfrac{D_2}{2} && \dfrac{E_2}{2} && F_2 \end{bmatrix}$
Points of intersection of the two equations will satisfy
$$ r^T Q_1 r = 0 $$
and
$$ r^T Q_2 r = 0 $$
Therefore, the points of intersection will satisfy
$$ r^T (\alpha Q_1 + \beta Q_2 ) r = 0 $$
Let $Q^* = \alpha Q_1 + \beta Q_2 $, we want to make the determinant of $Q^*$ zero for reasons that will become clear promptly. And to simplify the expression for the determinant, we'll consider first the case where $\alpha = 0$, if $\det(Q_2) = 0$, then we're done, otherwise, we can assume that $\alpha \ne 0$, and can therefore write
$Q^* = \alpha ( Q_1 + \dfrac{\beta}{\alpha} Q_2 ) $
Taking $\alpha = 1$ is as good as any non-zero value for $\alpha$, this simplifies $Q^*$ to
$Q^* = Q_1 + \beta Q_2 $
Since $Q_1$ and $Q_2$ are $3 \times 3$ symmetric matrices, then $\det(Q^*) $ will be a cubic polynomial in $\beta$, which must have at least one real root. Once we substitute this root into $Q^*$ we have our rank deficient matrix $Q^*$. A Rank deficient square matrix has at least one of its eigenvalues equal to zero, and this leaves the other two eigenvalues as both positive , both negative, or one positive and one negative. Diagonalizing $Q^*$ we obtain
$Q^* = R D R^T $
where $ D = \begin{bmatrix} D_{11} && 0 && 0 \\ 0 && D_{22} && 0 \\ 0 && 0 && 0 \end{bmatrix} $
Remember that our equation which we are trying to solve is
$$ r^T Q^* r = 0 $$
therefore, we want to solve
$$ r^T R D R^T r = 0 $$
Define $ w = R^T r $, i.e. $ r = R w $, then we end up with
$ w^T D w = 0 $
Recalling what $D$ is, this is
$ D_{11} w_1^2 + D_{22} w_2^2 = 0 $
Now if both $D_{11}$ and $D_{22}$ have the same sign, then the solution is $w_1 = w_2 = 0 $ while $w_3$ can be any real number. In this case, since $r = R w $ must have its third coordinate equal to $1$ , then this case corresponds to a single point.
We're not done, we have to verify that point we obtained satisfies
$ r^T Q_1 r = 0 $
If it does then this is our solution (a single point), otherwise there are no solutions.
Next, we'll consider the case where the eigenvalues have unequal signs, then we can choose that $D_{11} \gt 0 $ and $D_{22} \lt 0 $. In this case, the equation in vector $w$ becomes
$ D_{11} w_1^2 + D_{22} w_2^2 = 0 $
so that
$ w_2 = \pm \sqrt{ - \dfrac{ D_{11}}{D_{22}} } w_1 $
while $w_3 $ can be any real number. Thus in this case
$w = t (0, 0, 1) + s (1,\pm \sqrt{ - \dfrac{ D_{11}}{D_{22}} }, 0)$
That is $w$ lies in two possible planes spanned by the above two vectors for each of the two possible signs. Now
$ r = R w $
Again, a plane in the $r$ space for each choice of the sign in the above equation. Intersecting each plane (in $x, y, z$ of vector $r$ ) with the plane $z = 1$ (because we want the third coordinate of $r$ to be $1$ ), we obtain exactly two lines of intersection.
Each line is of the form
$ r = r_0 + \lambda d_0 , \hspace{10pt} \lambda \in \mathbb{R} $
For each of these lines, we now enforce
$ r^T Q_1 r = 0 $
This will ensure that $r^T Q_2 r = 0$ as well, because we know that $r$ satisfies $r^T (Q_1 + \beta Q_2) r = 0 $. This will give at most two solutions for the resulting quadratic equation in $\lambda$.
Thus we can have at most four solutions satisfying our system.
Here is a link to an Excel file that contains the source code for the function "solve-quadratic_system_in_two_unknowns3" that implements the above-mentioned algorithm to solve a system of two quadratic equations in two unknowns, in addition to many other related functions. These are included in this file as macros (VBA script). Click on the link, to open the online file, then click on "Editing" and choose "Open in Desktop App". This will open the file in your desktop Excel program. Click "View" then select "Macros".
It's not that hard to express $x$ in terms of $y$. But then if you don't know $y$, you might be stuck.
One possible approach to system of equation is elimination. That is, knowing that the equations in $x$ have a common solution, find a relation between the coefficients.
Say $a x^2 + b x + c = 0$ and $p x^2 + q x + r$ have a common solution $x$. This implies a relation between the coefficients $a$, $b$, $c$, $p$, $q$, $r$, see resultant. This can be calculated right away with some computer algebra system, say WA. From your system that would produce an equation in $y$ (of degree $4$ it turns out). Now, for each value of $y$, you get an $x$ ( in general one value only). This is all "in principle". Perhaps numerical methods might work better in concrete cases.
There is an alternate, elegant method of finding intersections of two conics. Homogenize and reduce to two quadratic equations $Q_i(x,y,z) = 0$, $i=1,2$. Now, determine $\lambda$ such that the conic $Q_1+ \lambda Q_2$ is singular ( this is an eigenvalue problem). For each of the three solutions $\lambda_i$, the conic $Q_1 + \lambda_i Q_2$ will be a product of two linear equations through the common points of the conics.
The gröbner basis of the intersection of two conics is generically of the form $\langle p_4(y),cx+q_3(y)\rangle$ (or $\langle p_4'(x),c'y+q_3'(x)\rangle$).
An example in M2:
R=QQ[x,y,MonomialOrder=>Lex]
I=ideal(x^2+29*x*y+23*y^2+19*x+17*y-13,x^2+2*x*y+3*y^2+5*x+7*y-11)
gens gb I -- matrix {{1507*y^4+3971*y^3-5041*y^2-7358*y-2012, 1318*x+40689*y^3+86119*y^2-179785*y-104956}}