since your function takes in all 3 coordinates, the only way to 'see' or plot the function is slicing the volume comprised by the domain used.
For such purpose the command slice comes really handy, have a look:
[x,y,z] = meshgrid(range,range,range);
v = x.*y.*z;
for k = range
hsp = surf(linspace(-N,N,numel(range)),linspace(-N,N,numel(range)), zeros(numel(range)) + k);
xd = hsp.XData;yd = hsp.YData;zd = hsp.ZData;
axis([-10 10 -10 10 -10 10])
note the handle h1 to pull the values on the slice plane, that may be useful to extract the slicing plane points with
the colour bar is also useful to visualise the actual value of the function.
if you find these lines useful would you please be so kind consider marking my answer as Accepted Answer?
To any other reader, if you find this answer of any help would you please click on the thumbs-up vote link,
thanks in advance for time and attention