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 $a\neq 0$ and $x\neq 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.