I'm trying solve a somewhat complicated system of PDEs using the Fourier transform. For the sake of completeness, the system is
\begin{align} \frac{\partial T}{\partial t}-K_x\frac{\partial^2 T}{\partial x^2}-K_z\frac{\partial^2 T}{\partial z^2} &= 0 \\ (1-2\nu) \left(\frac{\partial^2 u_x}{\partial x^2}+\frac{\partial^2 u_x}{\partial z^2} \right)+\frac{\partial^2 u_x}{\partial x^2}+\frac{\partial^2 u_z}{\partial x\partial z} &= 2(1+\nu)\alpha\frac{\partial T}{\partial x} \\ (1-2\nu) \left(\frac{\partial^2 u_z}{\partial x^2}+\frac{\partial^2 u_z}{\partial z^2} \right)+\frac{\partial^2 u_x}{\partial x\partial z}+\frac{\partial^2 u_z}{\partial z^2} &= 2(1+\nu)\alpha\frac{\partial T}{\partial z} \end{align}
where $u_{j} = u_{j}(t, x, z)$, and $T = T(t, x, z)$. The boundary conditions are
$$\lim_{z \to \infty} T(t, x, z) = 0, \qquad \frac{\partial T(t, x, 0)}{\partial z} = \frac{Q_0 \cos(q_0x)\delta(t)}{k_{z}}$$
while those for $u_x$ and $u_z$ are not relevant for this question.
Transforming to Fourier space ($x\rightarrow q$) and ($t\rightarrow \omega$) is pretty straightforward, and so is finding the solutions in Fourier space. I'm particularly interested in $u_z(t, x, 0)$, which in Fourier space is
$$\hat{u}_z(\omega,q,0)=\frac{A\eta(q)}{\beta(\beta+q)}$$
where $\eta(q)$ is a Fourier transform of $\cos(q_{0} x)$, which comes from the boundary condition for $T$. I packed the constants in $A$ and defined $\beta=\sqrt{\frac{q^{2} K_x + i \omega}{K_z}}$. Now I have to do inverse Fourier transforms ($q\rightarrow x$) and ($\omega\rightarrow t$). First, I do the inverse Fourier in time domain. I solve the integral using a Feynman trick, it's a bit lengthy so I won't go into detail, but the result is
$$\hat{u}_z(t,q,0)=B\eta(q)e^{-(K_x-K_z) q^{2} t} \text{erfc}(q\sqrt{K_zt})$$
Now, I can think of 3 approaches that in my mind should all lead to the same result.
$\eta(q)=\frac{1}{2}(\delta(q-q_0)+\delta(q+q_0))$, which after the inverse Fourier in $x$ gives $$u_z(t,x,z=0)=Ce^{-q_0^2(K_x-K_z)t}(\cos(q_0x)-i\sin(q_0x)\text{erf}(q_0\sqrt{K_zt}))$$
Let's replace $\cos(q_0x)$ in the boundary condition with $e^{iq_0x}$. In the end we will take a real part of the complex result. Then, $\eta(q)=\delta(q-q_0)$, which after taking only a real part leads to solution $$u_z(t,x,z=0)=Ce^{-(K_x-K_z)q_{0}^{2}t} \cos(q_0x) \text{erfc}(q_0\sqrt{K_zt})$$
Let's replace $\cos(q_0x)$ in the boundary condition with $e^{-iq_0x}$. In the end we will take a real part of the complex result. Then, $\eta(q)=\delta(q+q_0)$, which after taking only a real part leads to solution $$u_z(t,x,z=0)=-Ce^{-(K_x-K_z)q_{0}^{2} t}\cos(q_0x)(1+\text{erf}(q_0\sqrt{K_zt}))$$
So my question is, why does each of these approaches lead to a different solution and where exactly did I go wrong in the two that are incorrect. In the paper I'm studying, the authors don't explain the entire derivation, but the end result they state is the one from approach 2). It also does satisfy the initial equations.
Edit: Since there has been no answers or comments so far, I decided to share some of my findings. All three approaches lead to a solution that satisfies the initial system of PDEs. However, only approach 2 leads to a solution that also satisfies both boundary conditions. Approaches 1 and 3 do not satisfy $\lim\limits_{z \to \infty} u(t, x, z) = 0$, which was enforced early on by removing certain terms from the general solution. So somehow the approaches 1 and 3 reintroduced the terms that don't vanish at $z\rightarrow\infty$ and I would like to learn why that is so.