Perhaps translating mathematically will help clarify some things. For simplicity, I'll assume a 1D setup, with the EM fields and polarisation:
$$
\begin{align}
E_x(z,t) && B_y(z,t) && P_x(z,t)
\end{align}
$$
Maxwell's equations give:
$$
\begin{align}
\partial_zE_x &= -\partial_tB_y & -\partial_zB_y &= \mu_0\partial_tP_x+\frac{1}{c^2}\partial_t E_x
\end{align}
$$
Note that you can equivalently describe it with a single scalar quantity $A_x(z,t)$:
$$
\begin{align}
E_x &= -\partial_tA_x & B_y &= \partial_zA_x
\end{align}
$$
You'll need to add how to model matter:
$$
P_x = \epsilon_0\chi *E_x
$$
where $*$ is a temporal convolution product. I'll assume $\chi$ to be a Dirac for simplicity (the convolution is therefore a simple multiplication by a constant), but in general it is causal with a larger support.
Your plane of electrons can be modelled by taking:
$$
\chi = \chi_0\delta(z)
$$
so that the polarisation is localised at $z=0$. Now, adding an EM wave:
$$
\begin{align}
E_x^0 &= \hat E_xe^{i\omega(z/v-t)} & B_x^0 &= \hat B_xe^{i\omega(z/v-t)} \\
v &= \pm c & \hat B_y &= \frac{1}{v}\hat E_x
\end{align}
$$
From now on, the time dependence will always be $e^{-i\omega t}$ by linearity, so I will drop this factor.
The wave induces a polarisation:
$$
P_x = \epsilon_0\chi_0\hat E_x\delta(z)
$$
which in turn radiates waves (taking the retarded solution by causality):
$$
\begin{align}
A_x &= \int \frac{-i\omega}{(kc)^2-\omega^2}\chi_0\hat E_xe^{ikz}\frac{dk}{2\pi} \\
&= \frac{\chi_0}{2c}\hat E_xe^{i\omega|z|/c}
\end{align}
$$
which gives the fields:
$$
\begin{align}
E_x &= \frac{\chi_0i\omega}{2c}\hat E_xe^{i\omega|z|/c} \\
B_y &= \frac{\chi_0i\omega \text{ sgn}(z)}{2c^2}\hat E_xe^{i\omega|z|/c}
\end{align}
$$
As expected, you have two waves travelling away from the plane. Importantly, they both acquired a phase relative to the original wave.
The key fact is that in a medium, you have a continuum of polarisation planes. When a waves comes in, they radiate the previously calculated fields. However these scattered fields will in turn repolarise the medium and which will radiate light and so on ad infinitum. You have two possible approaches. Either solve the equations self consistently or sum up explicitly these scattered waves. The former is much easier and is done in all standard textbooks.
Taking back Maxwell's equations, you now assume $\chi$ to be independent on $z$. You still have:
$$
P_x = \epsilon_0\chi E_x
$$
but $E_x$ is the total EM fields which include the radiated fields due to the oscillation of $P_x$. You should therefore view this as a self-consistent equation. The usual method gives the dispersion relation with the index of refraction:
$$
n = \sqrt{1+\chi}
$$
This nonlinear dependence in $\chi$ translates the complicated resummation that has been bypassed.
If you truly want to understand how the individual phases conspire to a slow down the wave, you'll have to do the second method. The idea is to write the total field:
$$
\begin{align}
E_x &= \sum_nE_x^n & B_y &= \sum_nB_y^n
\end{align}
$$
with $E_x^0,B_y^0$ is your original plane wave and calculate by induction $E_x^{n+1}$ as the radiated waves by the polarisations induced by $E_x^n$. Using the previous formula and integrating over all the planes, you can write down the recurrence relation:
$$
\begin{align}
E_x^{(n+1)}(z) &= \chi\int\frac{i\omega}{2c}E_x^n(\zeta)e^{i\omega|z-\zeta|/c}d\zeta \\
B_y^{(n+1)}(z) &= \chi\int\frac{i\omega\text{ sgn}(z-\zeta)}{2c^2}E_x^n(\zeta)e^{i\omega|z-\zeta|/c}d\zeta
\end{align}
$$
As you can see, this is effectively a Taylor series in $\chi$. As Puk commented, it is important to take into account both the forward radiating waves and backward ones. Since we already know that the total field is:
$$
\begin{align}
E_x &= e^{i\frac{\omega z}{v}\sqrt{1+\chi}} & B_y &= \frac{1}{v}e^{i\frac{\omega z}{v}\sqrt{1+\chi}}
\end{align}
$$
Actually, using the previous remark, you can actually calculate the $E_x^n,B_y^n$ by Taylor expanding the formula in $\chi$. You'll therefore see how the phases interfere to give you a slowing down.
It is a nice complementary picture to the usual method, even if it's hard to do any explicit calculations with it. For example, one insight is that if $\chi$ is real and causal, the impulse response cannot travel faster than the speed of light, since it is true at every term of the perturbative expansion. You can check this explicitly with the solved expression, which is not completely obvious.
Caveat: the previous argument of causality only holds if the time response of the susceptibility $\chi(t)$ is regular enough at $t\to 0^+$, or equivalently that in the frequency domain $\hat \chi(\omega)$ decays for $\omega\to+\infty$. Indeed, say for a constant $\hat \chi \in (-1,0)$, the impulse response travels faster than the speed of light even if all the terms in the perturbative series travel slower than it. This is because the convergence is tricky to define in the time domain as it is essentially a Taylor expansion of a Dirac delta. The core of the argument is that formally:
$$
\delta(x+h) = \sum_{n=0}^\infty \frac{h^n}{n!}\delta^{(n)}(x)
$$
so even the distribution is supported on $\{-h\}$, you can write it as a "sum" of distributions supported on $\{0\}$. When trying to model materials more accurately, you do get a more regular susceptibility, so the impulse response will travel slower than the speed of light for realistic materials.
Hope this helps.
