>>9264161
First let us look at the symmetries of [math]\mathcal{L}[/math]. Suppose [math]T_\vec{a}[/math] and [math]T_\beta[/math] are representations of the space and time translation operators such that [math](T_\vec{a}f)(\vec{x},t) = f(\vec{x} - \vec{a},t)[/math] and [math](T_\beta f)(\vec{x},t) = f(\vec{x},t-\beta)[/math]. If [math][T_\vec{a},\mathcal{L}] = [T_\beta,\mathcal{L}] = 0[/math] then it follows that [math]G[/math] is also translationally invariant, and therefore [math]G(\vec{x},\vec{\xi},t,\tau) = G(\vec{x}-\vec{\xi},t-\tau)
[/math]. So you can WLOG solve the PDE [math]\mathcal{L}G(\vec{x},t) = \delta(\vec{x},t)[/math] to get the Green function, namely the Green function at the origin and [math]\tau=0[/math] determines the full Green function.
If the above is satisfied then you can perform a Fourier transformation. Given [math]G \in \mathcal{S}(\mathbb{R}^n)[/math] is a Riemann square-integrable distribution in space (i.e. [math]G\rightarrow 0[/math] as [math]\vec{x} \rightarrow \pm \infty[/math]), you can Fourier transform in the [math]\vec{x}[/math] coordinates and obtain an ODE in [math]t[/math] which you can solve depending on initial conditions, then inverse Fourier transform back into real space. If [math]G\in \mathcal{S}(\mathbb{M})[/math] is a square-integrable distribution in both the space and time coordinates then you can let [math]x = (t,\vec{x})[/math] and Fourier transform in [math]x[/math]. This gives you an algebraic equation that you can just inverse Fourier transform back directly.