Here is a heuristic algorithm which recognizes some (not all) surjective polynomials (this worked for me in practice).
The main idea is to try to find invertible polynomial map
$$ f, f_2 \ldots f_n \; : \mathbb{Q}^n \to \mathbb{Q}^n$$
Select bound $d$ for the degree of $f_2 \ldots f_n$
and make the coefficient of $f_i$ new variables $c_i$.
Compute the determinant $D$ of the jacobian matrix of $ f, f_2 \ldots f_n$
and try to solve symbolically for $c_i$, $D=1$.
If this succeeds, the jacobian conjecture implies the inverse
map is polynomial and solving the inverse map gives you solutions
as a side effect.
Added clarification answering Stefan's question
The range of $f$. It is $\mathbb{Q}$ as are the ranges of $f_i$.
In the example the given $f(x,y)$ is polynomial in x,y as is $f_2$.
In the example $A,B \in \mathbb{Q}$.
$f_i$ are auxiliary polynomials which are used by the jacobian conjecture
The coefficients of $f_i$. To construct the polynomials $f_i$,
for each $f_i$ generate all monomials in $x_i$ up to the chosen
degree $d$. $c_j$ are variables which are coefficients of each monomial in $x_i$, e.g. $c_{13} x_2 x_3$. So $f_i=\sum c_k \prod x_j$.
The determinant $D$ must be constant $\forall x_i$, so all coefficients
of $x_i$ except the constant must be $0$ and the constant coeff.
must be nonzero. In the given example, the solution allows some coefficients like $c_3$ to take any value.
About $c3 x$. This was copied from CAS and means $c_3 x^3$.
Example.
Take $f(x,y)={x}^{3}+3\,{x}^{2}y+3\,x{y}^{2}+{y}^{3}+3\,{x}^{2}+6\,xy+3\,{y}^{2}+2
\,x+3\,y$
Solving $D=1$ symbolicall gives
$$ f_2(x,y)= {\it c1}\,x+{\it c3}\,{x}^{3}+{\it c2}\,{x}^{2}+{\it c4}\,{x}^{4}+{
\it c20}\,{y}^{4}+{\it c15}\,{y}^{3}+{\it c10}\,{y}^{2}+{\it c5}\,y+{
\it c25}+{\it c24}\,{x}^{4}{y}^{4}+{\it c19}\,{x}^{4}{y}^{3}+{\it c23}
\,{x}^{3}{y}^{4}+{\it c14}\,{x}^{4}{y}^{2}+{\it c18}\,{x}^{3}{y}^{3}+{
\it c22}\,{x}^{2}{y}^{4}+{\it c9}\,{x}^{4}y+{\it c13}\,{x}^{3}{y}^{2}+
{\it c17}\,{x}^{2}{y}^{3}+{\it c21}\,x{y}^{4}+{\it c8}\,{x}^{3}y+{\it
c12}\,{x}^{2}{y}^{2}+{\it c16}\,x{y}^{3}+{\it c7}\,{x}^{2}y+{\it c11}
\,x{y}^{2}+{\it c6}\,xy$$
The inverse map of $f = A, f_2 = B$ is
$$
x=3\,{\it c3}\,A+3\,{\it c25}-A+{{\it c3}}^{3}{A}^{3}+{{\it c25}}^{3
}-{B}^{3}+6\,{\it c3}\,A{\it c25}-6\,{\it c3}\,AB-6\,{\it c25}\,B-6\,{
\it c3}\,A{\it c25}\,B-3\,B+3\,{{\it c3}}^{2}{A}^{2}+3\,{{\it c25}}^{2
}+3\,{B}^{2}+3\,{{\it c3}}^{2}{A}^{2}{\it c25}-3\,{{\it c3}}^{2}{A}^{2
}B+3\,{\it c3}\,A{{\it c25}}^{2}+3\,{\it c3}\,A{B}^{2}-3\,{{\it c25}}^
{2}B+3\,{\it c25}\,{B}^{2}$$
$$y=A-{{\it c3}}^{3}{A}^{3}-{{\it c25}}^{3}+{
B}^{3}-6\,{\it c3}\,A{\it c25}+6\,{\it c3}\,AB+6\,{\it c25}\,B+6\,{
\it c3}\,A{\it c25}\,B-2\,{\it c3}\,A+2\,B-2\,{\it c25}-3\,{{\it c3}}^
{2}{A}^{2}-3\,{{\it c25}}^{2}-3\,{B}^{2}-3\,{{\it c3}}^{2}{A}^{2}{\it
c25}+3\,{{\it c3}}^{2}{A}^{2}B-3\,{\it c3}\,A{{\it c25}}^{2}-3\,{\it
c3}\,A{B}^{2}+3\,{{\it c25}}^{2}B-3\,{\it c25}\,{B}^{2}
$$
This approach fails for $f = x y$ (modulo errors) and
succeeds for the Cantor pairing.
If you have specific examples, let me know to test my implementation.