The way it is implemented above, the value at 0 frequency is just the sum of your signal from -1 seconds to 1 second (you have about 5000 samples equal to 1 and about 5000 samples equal to 0)... so the expected amplitude should be about 5000.
A couple of suggestions that won't change much:
y = zeros(size(x));
y(x(i) >= -0.5 & x(i) <= 0.5) = 1;
and
Y = fft(y)/Fs; % Normalize FFT output
% to take into account that Fs is 5000, not 1
and
abs(Y) is equivalent to sqrt(Y.*conj(Y))
EDIT May 09, 2012
To plot something that looks more like a sinc, change your line:
plot(F,(sqrt(Y.*conj(Y))))
to
plot(F,real(Y))
and narrow your x-limits to +/- 6 (i.e., "xlim([-6 6])")
To confirm it is what you expect, create the sinc-like function you expect. From the "DFT pairs" link below, this would look something sort of like:
Z = zeros(size(Y));
Z(5001)=1;
Z(5002:5051) = -1*sin(pi*10000*(1:50)/20001)./(10000*sin(pi*(1:50)/20001));
Not sure why I needed the -1, but the positive frequency amplitudes seem to line up pretty well with what the expected shape is when you plot:
plot(F,real(Y),F,Z);
Just make sure you normalize the amplitudes in "Y" by Fs as I indicate above!
Also, to get an even time series change:
x=(-1:1/Fs:1)';
to
x=(-1:1/Fs:1-1/Fs)';
