Computing and Ploting Fourier transform

How can i use the fft command to transform a rectangular pulse to sinc function and plot the sinc function, i'm using a very traditional way to compute the fourier transform and plot it, but this way is too slow, when i use the fft command and try to plot, the ploting of the magnitude is only the pins or a spike, i need the plot to be a sinc function as the picture that i attached here, also i'm using the same method to compute the inverse fourier transform so i also need to use the ifft command, can anyone help?

if true
p = heaviside(t+1)-heaviside(t-1);
  k = 0;
F=zeros(1,8001);
for f= -4:0.001:4
  k=k+1;
F(k) = trapz(t,p.*exp(-1i*2*pi*f*t));
end
f = -4:0.001:4;
figure(1);
plot(f,F);
grid
end

Answer 1

[Rick Rosson]

Fs = 8;
dt = 1/Fs;
N = 8000;
t = dt*(-N/2:N/2-1)';

dF = Fs/N;
f = -Fs/2:dF:Fs/2-dF;

p = heaviside(t+1)-heaviside(t-1);

P = dt*fftshift(fft(ifftshift(p)));

figure;
plot(f,real(P));
grid on;

[omarahmed.m]

Thank you for your response, Could you just specify this command "P = dt*fftshift(fft(ifftshift(p)));", why it consist of fftshift of fft of ifftshift, and why multiplying by "dt". thank you very much.