Analytical solution for seismic wavefield in MATLAB
So, I know this question was asked a long time ago. But in the spirit of Stack Exchange, I will post my answer for future users perhaps.

The Greens functions (i.e., analytical solutions)

The first thing to be aware of are the various solutions to the wave equations change for the (1) source type and (2) the dimension.

The source type

Assuming acoustic waves in a homogeneous domain, that means you are either solving:
$$ \frac{\partial^2 p(x,t)}{\partial t^2} = c^2(x)\left(\nabla^2 p(x,t) +f(x,t) \right)  , \tag{1}$$
(where $p$ is pressure, $c$ the compressional velocity and $f$ the source function), or you are solving:
$$ \left\{ \begin{array}{rl} \frac{\partial p(x,t)}{\partial t} &amp; =  -c^2(x)\bigg( \nabla\cdot v(x,t) + g(x,t) \bigg), \\ \frac{\partial v(x,t)}{\partial t} &amp; = -\nabla p(x,t). \end{array} \right. \tag{2} $$
The difference between (1) and (2) is a time differentiation:
$$ f(x,t) = \frac{\partial}{\partial t}g(x,t). $$
This is an important difference in source functions!

The Greens function equation

It is standard to restate equation (1) in the following form:
$$ \left( \frac{1}{c^2(x)}\frac{\partial^2}{\partial t^2} - \nabla^2 \right)p(x,t) = f(x,t). \tag{3a}$$
Transforming this equation to the frequency domain, for which $\partial^2/\partial t^2 \to -\omega^2$, and where $\omega^2/c^2=k^2$ we obtain:
$$ \left( -k^2 - \nabla^2 \right)p(x,\omega) = f(x,\omega), \tag{3b}$$
or
$$ \left( k^2 + \nabla^2 \right)p(x,\omega) = -f(x,\omega), \tag{3c}$$
We can consult, e.g., https://www3.nd.edu/~atassi/Teaching/AME%2060633/Notes/greens.pdf to find the Greens functions (for which $f(x,\omega)=\delta(x-x_s)$ with $x_s$ the source location).

The Greens functions

I'll just copy the results. For equations of type (1) we have the Green's functions $G$:
$$
\begin{array}{lrl}
\text{1-D:} &amp; G &amp; = −\frac{i}{2k}e^{−ik\lvert x−x_s \rvert}, \\
\text{2-D:} &amp; G &amp; = −\frac{i}{4}H_0^{(2)}(k \lvert x−x_s \rvert),\\
\text{3-D:} &amp; G &amp; = \frac{1}{4\pi}\frac{e^{−ik\lvert x−x_s \rvert} }{\lvert x−x_s \rvert}.
\end{array}
$$
Indeed, $H_0^{(2)}$ is the Hankel function, in MATLAB we simply write that with besselh(0,2, k*r).

The MATLAB implentation

I made the following changes w.r.t. your code:


The frequency vector can be defined from $f_1\in(0,\frac{1}{\Delta t})$ rather than how you did it, which is $f_2\in(-\frac{1}{2\Delta t},\frac{1}{2\Delta t})$. In this way, the Hankel function only receives positive frequencies. Moreover, it is exactly the order in which the FFT produces its output. I thus also removed the fftshift operations.
I changed the Green's function to the 2-D case, i.e., I changed const=-1i*pi to const=-1i/4. 
I changed the nfft, to just using the FFT of the size of the wavelet or time vector. This made it easier to obtain results of similar amplitudes... Also, I don't divide the source function by nfft, to make the results relatively independent of the chosen dt.


Thus, the final code, having made only small changes compared to yours, becomes:

clc;    clear all;

%% Parameter Setting
vel=2000;   % velocity of medium
r=1000;     % distance at which wavefield will be observed
% the final wavelet should appear at time, t=r/vel=0.5sec

T=1;        % Total time
dt=.0025;   % time step
f0=25;      % central freq

%% Create Source Wavelet
N = round(T/dt);        % No of samples
t0 = 5/(sqrt(2)*pi*f0);     % Zero time shift
t = dt*(0:N-1);             % time vector
tau = t-t0+1*dt;     % time vector shifted by t0
src = (1 -2* tau.*tau * f0^2 * pi^2).*exp(-tau.^2 * pi^2 * f0^2);  % source wavelet

%%FT of signal
fmax=1/dt; % max frequency 
freq= linspace(0,fmax,N);    % Freq vector:   0 &lt;= f &lt;= Fmax
sw=fft(src);      % fourier transform of wavelet, S(w)

%% 2-D Greens Function 
w=2*pi*freq;
const = -1i/4;
H02 = besselh(0,2,w/vel*r);
H02(1)=0;
GF = const*H02;

%% Calculating the Wavefield  
WF= ifft(GF.*sw);

%% Plot
% Wavelet
subplot(2,2,1);     
plot(t,src);     title('\bf{Source wavelet}'); 
xlabel('Time(sec)');    ylabel('Amplitude ')

% Frequencies
subplot(2,2,2);   
plot(freq,abs(sw));    title('\bf{Source freq spec}');  
xlabel('Frequencies(Hz)');  ylabel('Amplitude ')
xlim([0 fmax/2])

% Gneens Fkt
subplot(2,2,3); 
plot(freq,abs(GF));   title('\bf{Green Function}'); 
xlabel('Frequencies(Hz)'); ylabel('Amplitude ')
xlim([0 fmax/2])

% Wafefield
subplot(2,2,4);      
hold on
plot(t,real(WF));    title('\bf{Trace}');   
xlabel('Time(sec)');  ylabel('Amplitude')


Provided that the signal is sampled to the Nyquist frequency, this should give good results: