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} & = -c^2(x)\bigg( \nabla\cdot v(x,t) + g(x,t) \bigg), \ \frac{\partial v(x,t)}{\partial t} & = -\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 <PERSON> 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 <PERSON> functions
I'll just copy the results. For equations of type (1) we have the <PERSON>'s functions $G$: $$ \begin{array}{lrl} \text{1-D:} & G & = −\frac{i}{2k}e^{−ik\lvert x−x_s \rvert}, \ \text{2-D:} & G & = −\frac{i}{4}H_0^{(2)}(k \lvert x−x_s \rvert),\ \text{3-D:} & G & = \frac{1}{4\pi}\frac{e^{−ik\lvert x−x_s \rvert} }{\lvert x−x_s \rvert}. 955
\end{array} $$ Indeed, $H_0^{(2)}$ is the <PERSON> 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:
1. 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 <PERSON> 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.
2. I changed the Green's function to the 2-D case, i.e., I changed const=-1i*pi to const=-1i/4.
3. 768