[Edited on 2013-10-02: see below for answer to the edited version of the question.]
In these sorts of problems it is always better if we know something about the system of equations. At the very least we should have some idea on how continuous the maps are (if they are not continuous then they are not computable and the whole thing is just unreasonable). Next, it is better to try to solve the system only in a bounded area (you can iteratively increase the area if needed).
Let us consider a simple example. Suppose we want to solve $\phi(\vec{x}) = 0$ on a unit cube in $\mathbb{R}^n$ up to precision $\epsilon$ and we know that $\phi$ is Lipschitz with constant $L$. Then we can sample $\phi$ on a grid of points such that no two are more than $\epsilon/L$ apart. We compute the sample values with precision $\epsilon/5$. Either we will find an approximate value $\vec{x}_0$ such that $|f(\vec{x}_0)| < 3\epsilon/5$, in which case we have an approximate solution, or all samples will have absolute values above $2\epsilon/5$ and there is no solution. This is going to be a very dumb and slow method (it is easy to write down Lipschitz functions with a huge $L$), but it is not easy to improve on it.
It really helps if $\phi$ is nice so that sensible numerical methods can be used. But in general the problem "find an $\epsilon$-zero or verify that all values are greater than $\epsilon/2$" is computable.
Answer to the edit in the question: If add the equation $\sin (\pi x_i) = 0$ to force $x_i$ to be integral, then the complexity is still quite horrible. We know that if there is a solution then it must be integral, so this might help us to quickly eliminate any solution that is obviously not integral, but it seems to be that in general complexity remains high.
As an example, consider the above situation of solving $\phi(\vec{x}) = 0$ on the unit cube in $\mathbb{R}^n$. The unit cube has $2^n$ points with integral coordinates. There does not seem to be a general method that will quickly eliminate those points, so naively we will expect complexity to be around $2^n \cdot t(\epsilon)$ where $t(\epsilon)$ is the complexity of evaluating $\phi$ up to precision $\epsilon$.
If we try to solve $\phi(\vec{x}) = 0$ on an unbounded region then the problem is not even computable, because (as you notice yourself) you can encode Diophantine equations in it.
So, if you have a cunning plan is to solve Diophantine equations, you should rethink the situation.