One way is to eliminate one of the variables (say, $y$) from the system, giving a quartic equation in $x$: namely, $$(31-x^2)^2 + x = 41 .$$ So solving the problem amounts to finding the roots of the quartic. [If you could guess one solution to the system, that of course helps in factoring the polynomial.]
–
SrivatsanDec 29 '11 at 16:23

As is already known to Lagrange, or even Euler, a system of two quadratic homogeneous equations define an elliptic curve on the plane. But your equation is not homogeneous. Therefore I might ask why did your friend ask this question? Thanks.
–
awllowerDec 29 '11 at 17:52

1

Then per chance you could try to solve the equations? $x^2+2y^2=1$ and $2x^2+y^2=1$ simultaneously. It is of course interesting.
–
awllowerDec 29 '11 at 17:56

This cubic polynomial is irreducible over the rationals or integers, which has three real solutions:

plot(x^3 + 5*x^2 - 37*x - 184,(x,-10,10))

For what it's worth, looking at the exact solutions, we can see that two are "variations" on the third, where the first term of each is multiplied by (the two nontrivial) third roots of unity, $e^{\frac{2{\pi}ik}{3}} = -\frac{1}{2}\pm\frac{i\sqrt{3}}{2},\;(k=1,2)$:

f.roots()[0] # latex(f.roots()[0])
# gives the first root with its multiplicity

A rational solution of the quartic equation has to be a divisor of $
920=2^{3}\times 5\times 23$. If we check $x=5$, we conclude that it is a
root, which we call $x_{0}$. Now we can easily factor the quartic
$$
x^{4}-62x^{2}+x+920=\left( x-5\right) \left( x^{3}+5x^{2}-37x-184\right) $$
remaining to solve the cubic equation
$$
x^{3}+bx^{2}+cx+d=0
$$
with coefficients $b=5,c=-37,d=-184$. To get the correspondent depressed cubic equation
$$
t^{3}+pt+q=0
$$
we need to make the change of variables $x=t-b/3$ thus finding the coefficients $p=-136/3$ and $q=-3053/27$. Now we can apply the Cardano's method. Since
the discriminant $q^{2}+4p^{3}/27<0$ all roots are real. One of
the roots is
$$
\begin{eqnarray*}
t_{1} &=&\left( -\frac{q}{2}+\frac{1}{2}\sqrt{q^{2}+\frac{4p^{3}}{27}}
\right) ^{1/3} +\left( -\frac{q}{2}-\frac{1}{2}\sqrt{q^{2}+\frac{4p^{3}}{27}}\right) ^{1/3}
\\
&\approx &7.742\qquad\text{(numerical evaluation in SWP)},
\end{eqnarray*}$$
which corresponds to $x_{1}=t_{1}-5/3\approx 6.075$. The radicals are chosen
in such a way that their product is

The remaining roots could also be found by factoring the cubic. Numerically we
have got $t_{2}\approx -4.487,t_{3}\approx -3.255$ and $x_{2}\approx
-6.154,x_{3}\approx -4.922$. The equation $y=31-x^{2}$ yields the values
$y_{0}=6, y_{1}\approx -5.910$, $y_{2}\approx -6.867$, $y_{3}\approx 6.776$ (Note: rounding error of $0.001$ in comparison with bgins's evaluation).

Rearranging gives $y^2 - x^2 = 10 + y-x,$ so that $(y-x) (y+x-1) = 10.$ You don't say you are looking for integer solutions, but if you are, then you are left with finitely many possibilities for the pair $(y-x,y+x-1),$ which can be individually solved by Gaussian elimination: for example, one possibility is $y-x = 2, y+x-1 =5,$ leading to $y = 4, x = 2.$ But this pair does not satisy the original equations. However, another possibility is $y-x = 1$ and $y+x-1 = 10,$ which leads to $y=6$ and $x=5,$ which does satisfy the original system. The other cases can be handled similarly.

Then this method is not sufficient, although the fact that $x=5$ is a solution of the quartic in Srivatsan's comment should be helpful.
–
Geoff RobinsonDec 29 '11 at 17:03

Notice that Srivatsan's quartic factors as $(x-5)(x^3 + 5x^2 -37x -184),$ so the remaining possibilities for $x$ that need to be checked (apart from $x=5$ which has already been dealt with) are the roots of the cubic factor.
–
Geoff RobinsonDec 29 '11 at 19:43

Note, if you graph $(31-x^2)^2$ again $41-x$ and $41$, you'll see that the constant line and the $41-x$ line are dominated by the large swing in the $(31-x^2)$, so the solutions to $(31-x^2)^2=41$ give you a good "guide" of where to look for the solutions of $(31-x^2)^2=41-x$. The solutions to the simpler equation are $\pm\sqrt{31\pm \sqrt {41}}$
–
Thomas AndrewsDec 29 '11 at 17:38

As a general algorithm for this kind of problem, you can compute a Gröbner basis of the ideal generated by your polynomials, using lexicographic ordering of monomials.

The result will be a "triangular system": first a polynomial in x, then a polynomial in xy... you can solve in x, then plug the values in the second polynomial and solve in y. (cf this paragraph in the WP article).

Solving Systems of Polynomial Equations (pdf) by Bernd Sturmfels gives an introduction to the topic. If you have access to the collective book Some Tapas of Computer Algebra (Springer), the two first chapters will give you a very complete and clear overview of the topic.