matlab - Fourier Transform of Morlet wavelet Function?

As you know the Morlet wavelet function is given by: $$\frac{1}{\sqrt{\pi f_b}}e^{\frac{-t^2}{f_b}}e^{j2\pi f_c}$$ The Fourier transform of this equation is: $e^{-\pi^2 f_b(f-f_c)^2}$ (is it right)?
First I attempted to plot the FFT of Morlet function by FFT function in Matlab, then I've plotted the Fourier transform function directly. I've expected to see same plots, but unfortunately they have different magnitude (why?). I publish the plot and my codes, do you have any idea why this happened? Is there any problem in my codes? Or any problem in the obtained Fourier transform? Or any problem in Matlab built-in FFT function?

I wrote this script to get those plots:

clear all;
fS = 500;
tStart= -4;
tStop= 4;
timeVector = linspace(tStart,tStop, (tStop-tStart)*fS );
fC = 2;
fB=2;
timeMask = zeros(1,length(timeVector));
timeMask((timeVector >= -fB/2) & (timeVector <= fB/2)) = 1;

psiWavelet = ((pi*fB))^(-0.5).*...                   % Morlet function
exp(2*1i*pi*fC.*timeVector).*exp(-timeVector.^2/fB); % Morlet function
%      FFT plot by matlab bulit-in FFT function
     Nfft =10*2^nextpow2(length(timeVector));
    FFT =fftshift(abs(fft(psiWavelet,Nfft)));
    freqs=[0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs=fftshift(freqs);
    figure(2);
    subplot(1,2,1)

    plot(freqs, FFT); 
    xlim([-1  5]);
    xlabel('Frequency / Hz');
    title (sprintf('Fourier Transform'));
%     FFT plot by its direct fourier transfrom function
    f_psi=exp(-(pi^2*fB)*(freqs-fC).^2);
    subplot(1,2,2)
    plot(freqs,f_psi)
    xlim([-1  5]);

Answer

If you approximate the Fourier transform

$$X(f)=\mathcal F(x)(f)=\int_{-\infty}^\infty x(t)\,e^{-2\pi j\,ft}\,dt$$

by the discrete Fourier transformation for by sampling on the time segment $[-T,T]$ as

$$X(f_n)\approx \sum_{k=-N}^{N-1} x(k\tau)\,e^{-2\pi j\,f_nk\tau}\,\tau=s[n]\,\tau$$

with $T=N\tau$, $f_n=n/(N\tau)=n/N*f_s=n/T$, $n=-N,...,N-1$, $s$ the result of FFT on the sampled $x$ sequence after shift, ...

then you have to multiply the result of the FFT method by $\tau=T/N=1/f_s$.