What I understood from your post is that you don't like the excess blur that happens at the downsampled image.
The theory of signal processing states that, when you downsample a digital signal by a factor of $M$ there's the potential of aliasing to happen if the signal is not bandlimited to $|w| < \pi/M$;
The spectral effect of the downsampling is such that the non-zero portions of the original signal spectrum above $ |w| > \pi/M$ can overlapp into fake frequency positions in the range $-\pi< w<\pi$ so that they will create visual artefacts that does not belong to the original signal.
In order to avoid such an overlapp (to avoid aliasing) the formal approach is to apply a lowpass filter with a cutoff frequency of $w_c = \pi/M$ and a gain of 1 at the passband.
The effect of this filter is to blur the image to the extend that when downsampled there won't be any significand overlap from the high frequency portions of the image. On the other hand this is eventually a blurring operation and it will throw away signal content (details of the image). So there will be a loss of quality (in terms of high frequency details) in your downsampled image.
However the problem in your code is that there is too much blurring, due to the Gaussian filter radius set to an inappropriate value. In theory what you need is an ideal filter with sharp transition at the cutoff frequency $w_c = \pi/M$ but the filters desiged in your code produces not only mild transition region but also a lower cutoff frequency than necessary.
So you can modify your lowpass filter radius untill you get an acceptable downsampled image.
Below is the modified code that shows spectrum of the Gaussian filter, you can use this to deduce if the filter cutoff is too much.
close all; clear all; clc;
I = double(imread('Cameraman.tif'))/255; % Black and Whit eimage in [0-255] range
figure,imshow(I);title('Original Image');
zoom_factor=2;
sigma=0.8*sqrt(zoom_factor-1); % I modified the radius here ???
[M,N,s]=size(I);
f1=-fix(M/2):ceil(M/2)-1;
f2=-fix(N/2):ceil(N/2)-1;
[f2,f1]=meshgrid(f2,f1);
Y=exp(-2*pi*pi*sigma*sigma*((f1/M).^2+(f2/N).^2)); % The BLUR filter
figure,imshow(abs(Y));title('Spectrum of the filter')
If=fft2(I);
If=fftshift(If);
for i=1:1
If(:,:,i)=If(:,:,i).*Y;
end
newI=ifftshift(If);
newI=real(ifft2(newI));
figure,imshow(newI);title('Gaussian Filtered Image')
newI=newI(1:zoom_factor:end,1:zoom_factor:end,:);
figure,imshow(newI);title('filtered and downsampled')
figure, imshow(newI);title('filtered and downsampled image')
figure, imshow(I(1:zoom_factor:end,1:zoom_factor:end));
title('downsampled without filtering')
Another solution is to use a highpass filter after the downsampling operation to improve the image by emphasizing the high frequency detail a bit.