This is hopeless without further assumptions, because no numerical procedure with a finite number of function evaluations can ever rule out a lack of injectivity in parts of the domain where the function hasn't been evaluated. In fact, even strengthening the hypothesis to assume that $f$ is smooth and taking $m=n=1$ isn't enough to make this possible.
Suppose that we have such a procedure, and start with a function $f$ that maps the interval [0,1] to [0,1] injectively. For example, $f(x)=x^{2}$ will work fine. Now suppose that we've evaluated the function at points
$x_{1}$, $x_{2}$, $\ldots$, $x_{n}$. Whatever these points are, we can always find an interval [a,b] that has been untouched. Now, let $g(x)=f(x)+w(x)$, where $w(x)$ is compactly supported on [a,b] but is a big enough "wiggle" to make $g(x)$ be no longer injective. Our procedure can't distinguish between $f(x)$, which is injective, and $g(x)$, which is not.
You might start to get somewhere by assuming smoothness and Lipschitz continuity of the derivatives. Assumptions like this are used in branch and bound algorithms for global optimization- there the goal is a to find a global minimum of a function that can only be evaluated numerically.