I am looking for an example (or a method to create my own examples) of diffeomorphism on a cube onto itself, similar to the wikipedia example on the page with the same name:
but in 3D. I'd like to have an example, but ideally I'd like to learn how to create examples myself. I understand the constrains that a diffeomorphism should have, and I think I'd be able to see if a given function is a diffeomorphism, but I dont know how to do the inverse.
(I guess that my problem is not a diffeomorphism)
I assume that you at least took a vector calculus and differential equations class (or classes). What you need is a smooth (at least continuously differentiable) vector field $V(x,y,z)$ on the cube $Q$, which vanishes on the boundary of the cube. For instance, you can (but do not have to) take $V$ as the gradient vector field of a smooth (at least twice continuously differentiable) function $f: Q\to {\mathbb R}$ whose 1st order partial derivatives vanish on the boundary of the cube:
$$
V=(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}).
$$
If you have hard time finding vector fields $V$, you can prescribe them numerically, on a sufficiently fine rectangular grid.
Once you have a vector field $V$, you consider the following ODE (ordinary differential equation; more precisely, a system of ODEs) $$ \frac{d}{dt} h(x,y,z,t)= V(x,y,z) $$ with the unknown function $h$ and the initial condition $h(x,y,z,0)=(x,y,z)$. Now, there is a theorem saying that this differential equation has a unique solution $h(x,y,z,t)$ for all $t\ge 0$, which defines a family of diffeomorphisms given, for each value of $t=t_0\ge 0$ (say, $t_0=1$), by the map $$ (x,y,z)\mapsto h(x,y,z, t_0). $$ Every such map will send the cube to itself, fix the boundary pointwise (i.e. send each boundary point to itself).
Solving ODEs is not easy but you should be able to find some software that does it for you numerically, which, I think, is all what you care about. For instance, both Mathematica and Matlab will do it. There is probably even some free software available.
You can also ask if every diffeomorphism of the cube to itself arises this way. This is not true but not far from being true. For instance, if you require that the diffeomorphism fixes some neighborhood of the boundary of the cube pointwise, then it will appear as a solution (for $t=1$) of a non-autonomous ODE very similar to the one I described above: $$ \frac{d}{dt} h(x,y,z,t)= V(x,y,z,t), $$ where $V$ is again a vector field on $Q$, but now it is also time-dependent ($V$ vanishes near the boundary of the cube). The initial condition is the same as above. The fact that all diffeomorphisms appear this way is a much-much harder theorem than the one I mentioned earlier. (For instance, this theorem is false once the dimension of the cube is $\ge 6$.)
Hope it helps.
There a second way (perhaps simpler) to construct lots of diffeomorphisms of the cube to itself, and with it you can probably more easily reconstruct the picture you posted in your question statement.
The basic idea is this:
If $f$ is a smooth function with domain the closed interval $(0,1)$, taking values in $\mathbb{R}$, such that $\lim_{x \to 0^+} f(x) = 0$ and $\lim_{x \to 1^-} f(x) = 1$, and $f' \geq \epsilon > 0$, then $f$ is a diffeomorphism of $(0,1)$ to itself.
The idea is that since $f$ is strictly monotonic, it is a bijection of $(0,1)$ to itself. And since $f' \neq 0$, it is everywhere a local diffeomorphism. Therefore $f$ is in fact a diffeomorphism.
Now, there is an easy way to guarantee that $f$ is a function with limits $\lim_{x\to 0^+} f(x) = 0$ and $\lim_{x\to 1^-} f(x) = 1$ and that $f' \geq \epsilon$: we can suppose that our function
$$ f(x) = x + g(x) $$
where $g(x)$ is any smooth function that vanishes at $x = 0$ and $x = 1$, and such that $|g'(x)| < 1 - \epsilon$.
This can be extended to the higher dimensional case as follows.
Let $\phi$ denote the desired diffeomorphism, and make the assumption that
$$ \phi(x_1, x_2, x_3) = (x_1,x_2,x_3) + (u,v,w) $$
where $u,v,w$ are real functions of $x_1, x_2, x_3$. We require that the three functions are smooth, that they vanish on the boundary of the cube, and that the sum of the gradients
$$ |\nabla u|^2 + |\nabla v|^2 + |\nabla w|^2 < 1 - \epsilon $$
This will guarantee that $\phi$ is a smooth map whose Jacobian derivative is invertible, and also using what is discussed in the first part of this answer that $\phi$ is bijective onto the cube.
Example:
In the two dimensional case, you can let
$$ u(x,y) = \frac{1}{4\pi} \sin(2\pi y) \sin(\pi x) $$
and
$$ v(x,y) = - \frac{1}{4\pi}\sin(2\pi x) \sin(\pi y) $$
and you will get a picture that looks something like what you showed above in the question.
In the cubic case, then you can do something similar by setting, for example,
$$ u(x,y,z) = \frac{1}{7\pi} \sin(\pi x) \sin(2\pi y) \sin(3\pi z) $$
$$ v(x,y,z) = - \frac{1}{7\pi} \sin(3\pi x) \sin(\pi y) \sin(2\pi z) $$
and
$$ w(x,y,z) = \frac{1}{7\pi} \sin(2\pi x) \sin(3\pi y) \sin(\pi z) $$
The crucial thing (which was used to determine the coefficient $\frac{1}{7\pi}$, is that we want the sum of the squares of the gradients of $u$, $v$ and $w$ to be strictly less than 1. A naive estimate using that trigonometric functions are bounded in absolute value by $1$ shows that $$ |\nabla u|^2 \leq \frac{\pi^2 + (2\pi)^2 + (3\pi)^2}{(7\pi)^2} = \frac{2}{7} $$ and similarly for $|\nabla v|^2$ and $|\nabla w|^2$, so the condition is verified.
