If $\hat u(\xi,t)$ is the Fourier transform of $u$, then the IVP becomes
$$
\hat u_t(\xi,t)=-\xi^2 \hat u(\xi,t) +\hat h(\xi,t), \quad \hat u(\xi,0)=0.
$$
Which implies that
$$
\hat u(\xi,t)=\int_0^t \mathrm{e}^{-\xi^2(t-s)}\hat h(\xi,s)\,ds.
$$
The only thing that remains is to use the inverse Fourier transform,
\begin{align}
u(x,t)&=\frac{1}{2\pi} \int_{-\infty}^\infty\int_0^t \mathrm{e}^{2\pi x\xi i-\xi^2(t-s)}\hat h(\xi,s)\,ds\, d\xi=\frac{1}{2\pi}\int_0^t\int_{-\infty}^\infty\mathrm{e}^{2\pi x\xi i-\xi^2(t-s)}\hat h(\xi,s)\,d\xi\,ds \\ &=\frac{1}{2\pi}\int_0^t\int_{-\infty}^\infty\mathrm{e}^{2\pi x\xi i-\xi^2(t-s)}\int_{-\infty}^\infty \mathrm{e}^{-2\pi i\xi yi}h(y,s)\,dy\,d\xi\,ds\\ &=\frac{1}{2\pi}\int_0^t\int_{-\infty}^\infty h(y,s)\left(\int_{-\infty}^\infty
\mathrm{e}^{2\pi x(\xi-y) i-\xi^2(t-s)}\,d\xi\right)\,dy\,ds.
\end{align}
It remains to show that
\begin{align}
\frac{1}{2\pi}\int_{-\infty}^\infty\mathrm{e}^{2\pi x(\xi-y) i-\xi^2(t-s)}\,d\xi=
\frac{1}{2\sqrt{\pi(t-s)}}\mathrm{e}^{-\frac{(x-y)^2}{4(t-s)}}.
\end{align}
The fourier transorm is typically taken with respect to the spacial variable, that is, the fourier transform of $f(x)$ is given by
$$\hat{f}(\xi) = \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty}f(x)e^{-i\xi x}\, dx$$
For us, it should be easy to show, using things like integration by parts and other integral techniques, that
$\widehat{A_{tt}(z,t)} = \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty}A_{tt}(z,t)e^{-i\xi z}\, dz = \hat{A}_{tt}(\xi, t)$
and
$\widehat{A_{z}(z,t)} = i\xi \hat{A}(\xi,t)$
Using these, we can transform your equation into the following:
$$i\omega a(\omega,t) = \frac{\alpha}{2}a(\omega,t)-\beta_1 a_{t}(\omega,t)-i\frac{\beta_2}{2}a_{tt}(\omega,t)$$
Now, this is an ordinary differential equation in t, so we may solve for $a(\omega,t)$ by solving this equation. Then, once we have an expression for $a(\omega, t)$, we can use the inverse fourier transform to write down the solution to the PDE.
Best Answer
You want to take the 2d FT over $x$ and $y$ (integration limits are $-\infty\to\infty$ in all integrals below): $$\hat{g}(k_x,k_y) \equiv \iint g(x',y')e^{-ik_xx'-ik_yy'}\,{\rm d}x'\,{\rm d}y' .$$ and with this convention the inverse transform is given by $$g(x,y) \equiv \frac{1}{(2\pi)^2}\iint \hat{g}(k_x,k_y)e^{ik_xx + ik_yy}\,{\rm d}k_x\,{\rm d}k_y .$$ This transforms your equation $u_{xx} + u_{yy} + k^2u = f(x, y)$ to $$\hat{u}(k_x,k_y)(-k_x^2-k_y^2+k^2) = \iint f(x',y')e^{-ik_xx' - ik_yy'}\,{\rm d}x'\,{\rm d}y'$$ Now dividing by $(-k_x^2-k_y^2+k^2)$ and taking the inverse transform you end up with the answer in the book (up to a renaming of the integration labels $x'\to \zeta$, $y'\to \eta$ etc.). Just be careful to not reuse an integration label already in use (this is a very common mistake and is why I used $x',y'$ instead of $x,y$ in defining the transform above).
When deriving this solution I made the implicit assumption that the FT exists (in the sense of a Riemann integral). When we take $f = 0$ we get $u = 0$ and this is indeed a solution. However this is only one of the solutions (and the same goes with the solution above). For example $u(x,y) = \sin(kx)$ is also a solution. The most general solution is obtained by adding the solution of the homogeneous equation $u_{xx} + u_{yy} + k^2u = 0$. None of the solutions of this equation has a FT (e.g. $\int_{-\infty}^\infty \sin(kx)e^{-ik_xx}{\rm d}x$ does not converge) in the traditional sense. The FT can be extended to a larger class of functions which has a FT in what we call a distribution (for which Dirac delta is one example) and we could make the derivation above work to produce these solutions, but that is a longer story (what changes is that an equation like $k\hat{u}(k) = \hat{f}$ no longer implies $\hat{u} = \frac{\hat{f}}{k}$ but rather $\hat{u} = \frac{\hat{f}}{k} + C\delta(k)$). Anyway separation of variables is a better way of solving the homogeneous equation if needed.