Step 1: Definition of multivariate derivative
By definition, $dg(x)$ is the unique matrix which satisfies:
$$
\lim_{h\rightarrow 0}\frac{||g(x+h)-(g(x)+dg(x) h)||_2}{||h||_2}=0
$$
In our case, this simplifies to:
$$
\lim_{h\rightarrow 0}\frac{||g(h)-h||_2}{||h||_2}=0
$$
Step 2: Working with cubes
But since we are working with cubes, not Euclidean balls, we would perfer the "cube-norm" (aka the sup-norm): $||x||_\infty:=\max_{i\in[n]}|x_i|$.
It still works, using the following inequalities (easy to show):
$$
||x||_\infty\leq ||x||_2\leq \sqrt{n}||x||_\infty
$$
We can write:
$$
0=\lim_{h\rightarrow 0}\frac{||g(h)-h||_2}{||h||_2}\geq \lim_{h\rightarrow 0}\frac{||g(h)-h||_\infty}{\sqrt{n}||h||_\infty}\geq 0
$$
So, it follows that the limit is zero, also for the $||\cdot||_\infty$-norm.
Finally, this means that for any $\varepsilon>0$, there exists a $\delta>0$, such that for all $||h||_\infty<\delta$ (in other words, for all $h\in C_\delta$),
$$
\frac{||g(h)-h||_\infty}{||h||_\infty}<\epsilon
$$
Step 3: Triangular inequality
Using the triangular inequality, $||h||_\infty\leq ||g(h)-h||_\infty+||g(h)||_\infty$.
Therefore:
$$
\frac{||h||_\infty}{||h||_\infty}-\frac{||g(h)||_\infty}{||h||_\infty}<\epsilon
$$
Rearranging:
$$
||h||_\infty(1-\epsilon)< ||g(h)||_\infty
$$
Using instead the triangular inequality $||g(h)||_\infty\leq ||g(h)-h||_\infty+||h||_\infty$, we get:
$$
\frac{||g(h)||_\infty}{||h||_\infty}-\frac{||h||_\infty}{||h||_\infty}<\epsilon
$$
Which means:
$$
||g(h)||_\infty\leq ||h||_\infty(1+\epsilon)
$$
To summarize:
$$
||h||_\infty(1-\epsilon)\leq ||g(h)||_\infty\leq ||h||_\infty(1+\epsilon)
$$
This already implies that $g(C_\delta)\subset C_{\delta(1+\epsilon)}$.
However, it is not enough to conclude that $g(C_\delta)$ contains all of $g(C_{\delta(1-\epsilon)})$.
Step 3: Inverse function theorem
By the inverse function theorem, there is some open set $V$ containing $0$, such that $g^{-1}$ exists on $V$ and is continuous.
Henceforth we will assume $C_{\delta(1+\epsilon)}\subset V$ (if not, we make $\delta>0$ smaller until this is true).
Moreover, the inverse function theorem tells us that $dg^{-1}(0)=I^{-1}=I$.
Repeat steps 1-2, to obtain:
$$
||g^{-1}(y)-y||_\infty\leq \epsilon||y||_\infty
$$
This holds for all $||y||_\infty<\gamma$ for some small $\gamma$, we may take $\gamma=\delta$ without loss of generality.
Step 4: Completing the proof
Take any $x\in C_{\delta(1-\epsilon)}$. We show that $x\in g(C_{\delta})$.
Since $g$ is invertible on $C_{\delta}$, this is the same as requiring that $y\in C_{\delta}$ whenever $g^{-1}(y)\in C_{\delta(1-\epsilon)}$.
Finally:
$$
||y||_\infty \leq ||g^{-1}(y)||_\infty+||g^{-1}(y)-y||_\infty\leq \delta(1-\epsilon)+\epsilon||y||_\infty
$$
Rearranging, we have $||y||_\infty<\delta$. The proof is complete.