Hi Andi,
The issue here is that for a voltage pulse the fourier transform is done on the linear voltage function y(t) to produce Y(f), but the FWHM condition is on the square of y(t) and Y(f). I modified your code to produce a gaussian pulse y(t) as the square root of your y(t) [which in this context is actually y(t)^2]. Then keep the cos(wt) carrier, do the fft to get Y(E), and plot y^2 and Y^2 for FWHM. The height of the peaks don't matter for FWHM, so I made no attempt to figure out the constants in front.
For the Y(E)^2 plot, the FWHM of each peak is about .09 eV (it's a bit undersampled), and 20fs times that is 1.8, so it works.
I added figure numbers since I'm not a fan of new figures popping up every time you rerun the code.
t=-100:0.01:100;
A=0.1519;
T0=20; %FWHM of the pulse
w=2*pi*300/800; %frequency centered at 800nm
y=sqrt(exp(-4*log(2)*(t/T0).^2)).*cos(w.*t); %Gauss pulse in terms of FWHM expression
y2 = y.^2;
figure(1)
plot(t,y2);
Y=fftshift(fft(y));
Y2_mag=abs(Y.^2);
X=(-(length(Y_mag)-1)/2:(length(Y_mag)-1)/2)*1239.84/300/length(Y_mag)/(t(2)-t(1)); %frequency axis
figure(2)
plot(X,Y2_mag,'-o');
xlim([1.4 1.7])
ylim([0 3e6])
