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.