Xrfgtjtk

I have the following simultaneous equations:

$$a x^2 + (b+2ay)x - c_1 = 0$$

$$ay^2 + (b+2ax)y - c_2 = 0$$

Where I'd like to solve for $x$ and $y$.

Obviously $a,b,c_1,c_2$ are known constants. They have the following restrictions:

$$a in [0,1]$$
$$b > 0$$
$$c_1, c_2 > 0$$

My thinking is that the above equations should give me a bunch of possible solutions, and then all but one can (hopefully) be eliminated according to restrictions that this particular problem imposes (eg. imaginary solutions are not permitted).

The only approach that I can see is to first use the quadratic formula to obtain solutions for $x$ in terms of $y$, and then substitute this into the second equation to then solve for $y$, or vise versa. But this quickly leads to an intractable (at least for me) expression. In particular, I end up with

$$x = frac{-(b+2ay)pm sqrt{(b+2ay)^2 + 4a c_1}}{2a}$$

along with the equation

$$- ay^2 pm sqrt{(b+2ay)^2+4ac_1}y - c_2 = 0$$

The hope is to then solve the above to obtain a solution for $y$, which can then be used to obtain the corresponding solution for $x$. But how on earth can we deal with the square root? Is this problem even solvable without resorting to a numerical solution?

You can make life simpler avoiding radicals but you will still arrive to a nasty quartic equation.

Considering $$a x^2 + (b+2ay)x - c_1 = 0 tag 1$$
$$ay^2 + (b+2ax)y - c_2 = 0 tag 2$$
eliminate $y$ from $(1)$to get (assuming $aneq 0$ and $xneq 0$)
$$y=frac{-a x^2-b x+c_1}{2 a x}tag 3$$
Plug this result in $(2)$, simplify as much as you can to get by the end
$$3 a^2 x^4+4a bx^3+(b^2+4 a c_2-2 a c_1 ) x^2 -c_1^2=0 tag 4$$

Now, have fun !

I think that considering some numerical method (such as Newton) could be a good idea looking for the zero's of function
$$f(x)=3 a^2 x^4+4a bx^3+(b^2+4 a c_2-2 a c_1 ) x^2 -c_1^2$$
Since $f(0)=-c_1^2 <0$, you could catch a first root (so you know that there is at least a second real root); deflate to a cubic equation (using long division) and this would be easier to handle.

Edit

Because of $f(0)=-c_1^2 <0$, we can think about applying Newton method twice first. Taking advantage of $f'(0)=0$, the estimates of two roots can be obtaines using a series expansion around $x=0$. This will give as estimates
$$x_pm=pm sqrt{frac{c_1^2}{b^2+4 a c_2-2 a c_1 }}$$

Let us try using $a=3$, $b=7$, $c_1=31$, $c_2=50$. This would give $x_pm=pm frac{31}{sqrt{463}}approx pm 1.44$. Let us start Newton method for the zero's of function
$$f(x)=27 x^4+84 x^3+463 x^2-961$$
For the positive root, the iterates would be
$$left(
begin{array}{cc}
n & x_n \
0 & 1.440692178 \
1 & 1.272119191 \
2 & 1.254474954 \
3 & 1.254291856 \
4 & 1.254291837
end{array}
right)$$
For the negative root, the iterates would be
$$left(
begin{array}{cc}
n & x_n \
0 & -1.440692178 \
1 & -1.559623512 \
2 & -1.554550306 \
3 & -1.554540654
end{array}
right)$$ Now, making the division
$$frac{27 x^4+84 x^3+463 x^2-961 }{(x-1.254291837)(x+1.554540654)}=27 x^2+75.8933 x+492.859$$ which does not show any real root. So, we have our solutions for $x$ and back to $(3)$ the corresponding $y$'s.

Qualitatively, Bézout’s Theorem on plane curves says that there will be four intersection points, including (pairs of conjugate) complex points and points at infinity.

A glance at the quadratic terms ($ax^2+2ayx$ in the first equation, $ay^2+2ayx$ in the second) shows that there will be no intersection points at infinity. In fact, the curves are both hyperbolas, with entirely different asymptotes. Thus you expect two pairs of complex points (no real intersection), one such pair and two real points (two visible points of intersection), or four visible intersection points, always taking account of multiplicity of course.

In the event, thanks to Desmos, it looks as if there are always at least two visible points of intersection; certainly when you take $a=1$, $b=4$, $c_1=1$ and $c_2=frac16$, you get four honest real points of intersection.