I will give you an answer for the part where N is even. The rest can be easily adapted, I hope you can do this yourself :-) Or at least try it - if you've got problems, just come back to us.
I don't have MATLAB installed anymore, so I hope this is free of typo errors. But the idea should be clear:
%
N = 10;
% First create all possible starting coordinates of 2x2x2 cubes within the big cube:
coords = 1:2:(N-1); % could easily be adapted to other small-cube sizes like 3x3x3 if you want to...
% Now create all possible combinations of starting-coordinates in every direction (as it is a cube, the starting points in x, y, z directions are the same):
sets = {coords, coords, coords};
[x y z] = ndgrid(sets{:});
cartProd = [x(:) y(:) z(:)]; % taken from here: http://stackoverflow.com/a/4169488/701049 --> you could instead also use this function: https://www.mathworks.com/matlabcentral/fileexchange/10064-allcomb-varargin- which generates all possible combinations
% Now cartProd contains all possible start-points of small cubes as row-vectors. If you want, you can easily calculate the corresponding vectors of end-points by simply adding +1 to every entry which will effectively yield a small-cube size of 2. If you want to further modify it to support e.g. 3x3x3 cubes, simply add +2 to every dimension.
endPoints = cartProd+1;
% E.g.: the first small cube starts at [x,y,z] = cartProd(1,:) and ends at [x_e, y_e, z_e] = endPoints(1,:).
Have fun :-)
Hint: for the odd big cube -> Simply treat it as evenly-sized cube, e.g. treat a 9x9x9 cube as 10x10x10, take my algorithm from above, and then move the most outer small-cubes one step to the center. That means, take the small cubes with the biggest x, y or z coordinate and substract 1 in that direction. So that the starting coordinate for all small cubes with x_max=9 will be changed to x=8. Then the same for y_max=9 and z_max=9.