If you have the Image Processing Toolbox, I would make this problem easier on yourself by using imwarp and possibily imref3d. You can represent what you are trying to do as a geometric transformation of your input image with a scale transformation. It sounds like your intention is to choose a scale factor > 1 to increase the grid size of your image.
Using imref3d will help you define the world location of your input image from the metadata in your dicom image.
S = load('mri');
I = squeeze(S.D);
RI = imref3d(size(I),0.9375 0.9375 1.200);
scaleFactor = 2;
tformResize = affine3d([scaleFactor 0 0 0; 0 scaleFactor 0 0; 0 0 scaleFactor 0; 0 0 0 1]);
[Inew,Rnew] = imwarp(I,RI,tformResize);
Conceptually, it should be possible to do what you are trying to do with interp3. My guess is that interp3 is returning a single/double image outside the range [0 1] because you are probably just using simple cast single(I) or double(I) instead of rescaling and casting (im2double). Many MATLAB/IPT functions that require a dynamic range assume that single/double data is normalized to the range [0 1]. I bet you have not actually lost information if you are interpolating correctly. You either need to adjust the dynamic range used in visualizing your floating point data:
imshow(I(:,:,10),[]);
Move your output image back to its original datatype, for example uint8, or rescale your output image to the range [0 1]:
I = I./max(I(:));
