Ok, figured it out: ifftshift (or fftshift) is only for plotting, so before doing the fft (or ifft) the complex function has to be regenerated by ifftshift (or fftshift) - here is the code:
%-----------------------------------
w1 = linspace(0,1000,1024);
w1 = w1;
W = 10;
w0 = 100;
N = length(w1);
dw = mean(diff(w1));
Y1 = -exp(-((w1-w0)/W).^2);
%PLOT
subplot (1,3,1); %plot of the original frequency domain signal I started with
plot (w1,Y1)
xlim ([-200 200]);
axis square;
%IFFT
y1 = ifft(Y1);
y1 = ifftshift(y1); %this is only for plotting
t1 = (-N/2 : N/2 - 1)/dw;
%PLOT
subplot (1,3,2);
plot (t1, real(y1)) %plot of the time domain signal
xlim ([-100 100]);
axis square;
%FFT
y1 = ifftshift(y1); %this is crucial
Y2 = fft(y1);
Y2 = fftshift(Y2); %this is only for plotting
w2 = (-N/2 : N/2 - 1)*dw;
%PLOT
subplot (1,3,3);
plot (w2,real(Y2)) %plot of the real part of frequency domain signal
xlim ([-200 200]);
axis square;
%-----------------------------------