Correct scaling for FFT to mimic CT FT?

Hi, I've got myself confused with the FFT implementation of the analytical continuous time FT, and wondered if anyone could help.

For a square pulse of width tau seconds and amplitude Amp, the analytical Fourier Transform is G(jf) = Amp*tau*sinc(f*tau) = Amp*tau*sin(pi*f*tau)/pi*f*tau. However when I carry this out using the FFT (using the code below), I get a different amplitude of G(jf).

I have also attached an image depicting results. The pulse has width 4 (centered at t = 0), and amplitude 13.2. The analytical FT expression is G(jf)=13.2*4*sinc(f*tau) which has a maximum amplitude of 52.8, but my implementation using the FFT has a maximum amplitude of 52,813 (1000 times greater). I've tried various scaling factors related to the length of the original vector N etc, but still can't work it out - what should one scale by to give the correct data corresponding to the analytical FT? and why?

Kind regards, Olie.

Code below:

%%Time vector:
dt = 0.001;
t = -10:dt:10;

%%Function g (pulse):
g = zeros(length(t),1);
Amp = 13.2;
g(8001:12001) = Amp; % amplitude
tau = 4;             % width
figure(1)
subplot(3,1,1)
plot(t,g); grid on; xlabel t; ylabel g(t); title 'original function, g(t)'
string1 = ['amplitude, Amp = ' num2str(Amp)];
text(3,10,string1)

%%Actual FT:
freq = -3:0.0001:3;
Gact = Amp*tau*sinc(tau.*freq);
figure(1)
subplot(3,1,2)
plot(freq,abs(Gact)); grid on; xlabel 'freq / Hz'; ylabel |G(jf)|; title 'actual FT, G(jf) = Amp*tau*sinc(f*tau)';
string1 = ['amplitude = ' num2str(max(Gact))];
text(0.5,40,string1)

%%FT using normal FT code (discrete_fourier_transform):
% Number of points (N-point DFT)
N = length(g);

% Range of t
t = [0:dt:(N-1)*dt];

% DFT data
Gnormal = fft(g);
Gnormal = Gnormal(1:ceil(N/2+1),1);  % Truncates data to only one half (right half is mirror)

% Range of discrete frequencies (rad/sample)
omega_dis = [0:2*pi/N:ceil(N/2)*2*pi/N];

% Range of continuous physical frequencies (rad/s)
omega_cont = omega_dis/dt;

% Range of continous physical frequencies (Hz)
freq = omega_cont/(2*pi);

figure(1)
subplot(3,1,3)
plot(freq,abs(Gnormal)); grid on; xlabel 'freq / Hz'; ylabel |G(jf)|; title 'FFT, G(jf)';
string1 = ['amplitude = ' num2str(max(real(Gnormal)))];
text(0.5,4e4,string1)
axis([-3 3 -1000 7e4])