Finally, I've accomplished my goal by individually `surf`

-plotting each face of the transformed cube:

After looking at the docs of `surf`

I found `sphere(n)`

to produce an input for `surf`

which would render the unit sphere and looked into it: `surf(X,Y,Z)`

where all are matrices of the same size produces the surface with these points preserving the rectangular topology. This leads to the following code:

```
[e1x, e1y] = meshgrid(tx,ty);
[e2x, e2z] = meshgrid(tx,tz);
[e3y, e3z] = meshgrid(ty,tz);
% These are the three different grids needed for the cube faces
e1x = e1x(:);
(...)
% Make all of them columns
e1z = ones(size(e1x));
(...)
% The respective third dimensions will be constant 0 or 1
bottom = [ e1x e1y 0*e1z]';
top = [ e1x e1y 1*e1z]';
front = [ e2x 0*e2y e2z]';
back = [ e2x 1*e2y e2z]';
left = [0*e3x e3y e3z]';
right = [1*e3x e3y e3z]';
% Create the faces:
% surf(bottom(1,:), bottom(2,:), bottom(3,:)) with appropriate reshape
% would plot the bottom face of the unit cube
x = [bottom top front back left right];
Kstart = cumsum([1 size(bottom, 2) size(top, 2) size(front, 2) size(back, 2) size(left, 2)]);
Kend = [Kstart(2:end)-1 size(x,2)];
% Concatenate and save the start- and end-indices of each face
Ndim = [Ny Nx;Ny Nx;Nz Nx; Nz Nx; Nz Ny; Nz Ny];
% Dimensions for each face
```

And then after evaluation of `x`

into the variables `Px, Py`

and `Pz`

:

```
hold on
for k=1:length(Kstart)
Kcurr = Kstart(k):Kend(k);
% Index range of the points for the k-th face
surf(reshape(Px(Kcurr), Ndim(k,:)), ...
reshape(Py(Kcurr), Ndim(k,:)), ...
reshape(Pz(Kcurr), Ndim(k,:)), ...
k*ones(Ndim(k,:)), 'EdgeColor', 'none');
% surf each face assigning color k to face k for distinction
end
```

This produces the following image:

Note that the red line is where two faces meet inside the figure.