I have came recently across the following equation
$$\psi(x,t_2)=\int G(x,y)\psi(y,t_1)dy$$
I want to understand its interpretation. Here is what i understand, I see that this equation gives the form of the new wavefront (which i believe is a straight line) $\psi(x,t_2)$ at the new position $x$ and the new time $t_2$ given the integral of the previous small wavelets multiplied by the propagator $G$ over the entire y axis, that is given $\int G(x,y)\psi(y,t_1)dy$
What you are looking at is the representation of a linear operator as in integral (which is not possible for all of them).
More specifically, the operator you are looking at, is a time-translation invariant, time-evolution operator $U(t_2, t_1)$. This operator evolves your solution at time $t_1$ to your solution at $t_2$, in other words $$ \psi(x, t_2) = \big(U(t_2, t_1) \psi\big)(x, t_1) $$
If such an operator is reasonably nice, it can be represented as an integral (note that the word "operator" typically implies linearity): $$ \psi(x, t_2) = \int dx'\, G(x, x') \psi(x', t_1) $$
The label $y$ is just chosen randomly in your original equation (I have chosen $x'$), it's just a dummy integration variable, you can label it as you want). $x$ need not even be a number in this notation, it could be an element of Euclidean space $\mathbb{R}^n$ and the equations would all look the same. It refers to the same physical space as the variable $x$ of the resulting $\psi(x, t_2)$.
$G(\cdot, x)$ can be understood as the wave after a evolution time of $t_2 - t_1$ when the initial conditions where $\psi(x, t_1) = \delta(x)$. With this we can make your intuitive description of the structure of the solution explicit, here we write $U(t_2, t_1)$ as $U_x$ to make clear that it only acts on functions $x \mapsto f(x)$ in the free variable $x$: $$ \psi(x, t_2) = U_x \psi(x, t_1) = U_x \int dx' \psi(x', t_1) \delta(x - x') = \int dx' \psi(x', t_1) U_x \delta(x - x') = \int dx' \psi(x', t_1) G(x, x'). $$ So the solution $\psi(x, t_2)$ is the superposition of the time-evolution of the delta spikes weighted by the initial state. We could pull the operator into the integral here since the operator is linear (and we assume it and our initial data to be "reasonable nice") and $U_x$ does not act on $x'$ so the $\psi(x', t_1)$ are just coefficients.
Note, that this is not the Huygens principle as usually given – the solution denpends not only on the points on the wavefront but on the value of the solution in the entire support of $G$ at an earlier time. (The wave equation $\partial_t^2 \phi = c^2 \Delta \phi$ in odd dimensions, however, does permit a solution that can be interpreted as strictly following the Huygens principles).
In other words: In general, you don't construct wavefronts from wavefronts, but rather wave configurations from wave configurations. And in this sense it is a strongly generalized Huygens principle.
Note: Wikipedia has quite a nice discussion of more or less exactly these issue: https://en.wikipedia.org/wiki/Huygens%E2%80%93Fresnel_principle.
