The answer is 'no', because polynomial mappings with polynomial inverses preserve volumes up to a constant multiple.

To see why this property holds, suppose that $p:\mathbb{R}^d\to\mathbb{R}^d$ is a polynomial mapping with polynomial inverse $q:\mathbb{R}^d\to\mathbb{R}^d$. Then $p$ and $q$ extend to $\mathbb{C}^d$ as polynomial maps with polynomial inverses. This means that the Jacobian determinant of $p$ on $\mathbb{C}^d$ is a complex polynomial with no zeros and hence must be a (nonzero) constant.

Now, consider a diffeomorphism $f:\mathbb{R}^d\to\mathbb{R}^d$ that is *radial*, i.e., $f(x) = m(|x|^2)x$ for some smooth function $m>0$. One can easily choose $m$ in such a way that $m(4)=1/2$ and $m(9)=4/3$, so that $f$ maps the ball of radius $2$ about the origin diffeomorphically onto the ball of radius $1$ about the origin while it maps the ball of radius $3$ about the origin diffeomorphically onto the ball of radius $4$ about the origin.

Let $\epsilon>0$ be very small and suppose that $\|f-p\|_{\infty;U} <\epsilon$ for $U$ chosen to be some very large ball centered on the origin. Then $p$ maps the sphere of radius $2$ about the origin to within an $\epsilon$-neighborhood of the sphere of radius $1$, while it maps the sphere of radius $3$ about the origin to within an $\epsilon$-neighborhood of the sphere of radius $4$. It's easy to see from this that $p$ cannot have constant Jacobian determinant.

**Added remark:** The group $\mathrm{SDiff}(\mathbb{R}^d)$ consisting of volume-preserving diffeomorphisms of $\mathbb{R}^d$ is a 'Lie group' in Sophus Lie's original sense (i.e., a group of diffeomorphisms defined by the satisfaction of a system of differential equations; in this case, that the Jacobian determinant be equal to $1$).

The subgroup $\mathcal{SP}(\mathbb{R}^d)\subset \mathrm{SDiff}(\mathbb{R}^d)$ consisting of volume-preserving polynomial diffeomorphisms with polynomial inverses however, is *not* a 'Lie subgroup' in Lie's original sense when $d>1$, as it cannot be defined by the satisfaction of a system of differential equations: It contains all of the mappings of the form $p(x) = x + a\,(b{\cdot}x)^m$ where $a,b\in\mathbb{R}^d$ satisfy $a\cdot b = 0$ and $m>1$ is an integer (indeed, $p^{-1}(y) = y - a\,(b{\cdot}y)^m$), plus, it contains $\mathrm{SL}(d,\mathbb{R})$ and the subgroup consisting of the translations. Using this, it is easy to show that, for any $f\in\mathrm{SDiff}(\mathbb{R}^d)$ and for any integer $k$, there exists a $p\in \mathcal{SP}(\mathbb{R}^d)$ such that $f$ and $p$ have the same Taylor series at the origin up to and including order $k$. Thus, $\mathcal{SP}(\mathbb{R}^d)$ cannot be defined by a system of differential equations (in Lie's sense).

Using this Taylor approximation property, one can prove that $\mathcal{SP}(\mathbb{R}^d)$, like $\mathrm{SDiff}(\mathbb{R}^d)$, acts transitively on $n$-tuples of distinct points in $\mathbb{R}^d$ for any integer $n$. Whether one can prove that $\mathcal{SP}(\mathbb{R}^d)$ can 'uniformly approximate' $\mathrm{SDiff}(\mathbb{R}^d)$ on compact sets is an interesting question.

1more comment