Edit :
I tried to complete the Feynman calculation. Thank you for being kind to English, which is not my native language ! And sorry for possible miscalculation !
The problem is part of the Oseen extinction theorem.
You'll find usueful reference in https://en.wikipedia.org/wiki/Ewald%E2%80%93Oseen_extinction_theorem
In particular, the irremplace Born and Wolf. But it's complicated !
This is a simpler calculation. It is a macroscopic one. It is not perfect but I think it is interesting.
The field in a point is the sum of all the fields radiated by all the slices. But these radiated fields depend of the fied himself and so we are led to an integral equation.
In order not to be too long, I consider as known the field radiated by a plane of uniform polarization which oscillates sinusoidally.
A thin plate placed at $z=0$ and of length $dz$ is polarized : a volume $\text{d}\tau $ is equivalent to the dipole moment $\overrightarrow{\text{d}p}=\overrightarrow{P}\text{d}\tau =\overrightarrow{P}\text{d}z\text{d}S$. The polarization is supposed to vary sinusoidally in time and directed according to Ox: $\overrightarrow{P}=\overrightarrow{{{P}_{0}}}{{\operatorname{e}}^{j\omega t}}$with $\overrightarrow{{{P}_{0}}}={{P}_{0}}\overrightarrow{{{e}_{x}}}$
We work in complex representation but the complex amplitudes are not underlined to lighten the notations. All the vectors are along Ox.
The field radiated by this plate is ($k=\omega /c$):
$\text{d}{{E}_{x}}=-\frac{1}{2{{\varepsilon }_{0}}}(jk){{P}_{0}}dz{{\operatorname{e}}^{j\left( \omega t-kz \right)}}$ if $z>0$ and $\text{d}{{E}_{x}}=-\frac{1}{2{{\varepsilon }_{0}}}(jk){{P}_{0}}dz{{\operatorname{e}}^{j\left( \omega t+kz \right)}}$ if $z<0$
On can find the result in the Feynman's lectures : http://www.feynmanlectures.caltech.edu/I_30.html
I have simply adapted this result to a polarized medium.
Now consider a medium such that the polarization ${{P}_{x}}$ is related to the complex amplitude of the electric field by the relation ${{P}_{0}}={{\varepsilon }_{0}}\chi {{E}_{x}}$with the a priori complex quantity called susceptibility and dependent (a priori) of the frequency. The complex indices is linked to the susceptibility by ${{n}^{2}}=1+\chi $ . (This definition is easily adapted to a conductor with conductivity $\gamma $ )
The medium extends for z varying from 0 to infinity. For, $z<0$ we have vacuum.
The medium is "lit" by an incident plane wave polarised along Ox, ${{E}_{i}}={{E}_{0}}{{\operatorname{e}}^{j}}^{(\omega t-kz)}$ with $k=\omega /c$. This wave will polarize the medium as in the previous question.
Each slice of the medium will then radiate. The total field at a point z, is the sum of the incident field and the fields radiated by the various slices.
By changing the origin for each slice, we have :
${{E}_{z'>z}}(z)=-\frac{1}{2{{\varepsilon }_{0}}}jk\int\limits_{z}^{+\infty }{{{P}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t+k(z-z') \right)}}}$
${{E}_{0<z'<z}}(z)=-\frac{1}{2{{\varepsilon }_{0}}}jk\int\limits_{0}^{z }{{{P}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t-k(z-z') \right)}}}$
Since ${{P}_{x}}(x')={{\varepsilon }_{0}}\chi {{E}_{x}}(x')$ :
${{E}_{x\,}}_{z'>z}(z)=-\frac{1}{2}jk\chi \int\limits_{z}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t+k(z-z') \right)}}}$
${{E}_{x\,0<z'<z}}(z)=-\frac{1}{2}jk\int\limits_{0}^{z }{\chi {{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t-k(z-z') \right)}}}$
At $z>0$ , the amplitude of the total field satisfies the integral equation:
${{E}_{x}}(z,t)={{E}_{0}}{{\operatorname{e}}^{j(\omega t-kz)}}\underbrace{-\frac{1}{2}jk\chi \int\limits_{z}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t+k(z-z') \right)}}}}_{{{E}_{x\,}}_{z'>z}(z)}\underbrace{-\frac{1}{2}jk\chi \int\limits_{0}^{z}{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t-k(z-z') \right)}}}}_{{{E}_{x\,}}_{z'<z}(z)}$
Since the medium is linear, we have ${{E}_{x}}(z,t)={{E}_{x}}(z){{\operatorname{e}}^{j(\omega t)}}$
and the field obey the integral equation, valid if $z>0$ :
${{E}_{x}}(z)={{E}_{0}}{{\operatorname{e}}^{-jkz}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{jkz}}\int\limits_{z}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{-jkz'}}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{-jkz}}\int\limits_{0}^{z}{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{+jkz'}}}$
To solve this integral equation, we can try a solution of the form ${{E}_{x}}(z)=C{{\operatorname{e}}^{\left( -j\beta z \right)}}$with C and $\beta $ constants a priori complex, and with $\operatorname{Im}(\beta )>0$ to avoid divergence. (We could also differentiate the equation twice).
$\int\limits_{0}^{z}{\text{C}{{\operatorname{e}}^{-j\beta z'}}\text{d}z'{{\operatorname{e}}^{+jkz'}}}=C\frac{1}{j(k-\beta )}({{\operatorname{e}}^{+j(k-\beta )z}}-1)$
$\int\limits_{z}^{\infty }{\text{C}{{\operatorname{e}}^{-j\beta z'}}\text{d}z'{{\operatorname{e}}^{-jkz'}}}=-C\frac{1}{j(k+\beta )}({{\operatorname{e}}^{-j(k+\beta )\infty }}-{{\operatorname{e}}^{-j(k+\beta )z}})=C\frac{1}{j(k+\beta )}{{\operatorname{e}}^{-j(k+\beta )z}}$
If we replace in the original equation:
$C{{\operatorname{e}}^{-j\beta z}}={{E}_{0}}{{\operatorname{e}}^{-jkz}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{jkz}}C\frac{1}{j(k+\beta )}{{\operatorname{e}}^{-j(k+\beta )z}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{-jkz}}C\frac{1}{j(k-\beta )}({{\operatorname{e}}^{+j(k-\beta )z}}-1)$
$C{{\operatorname{e}}^{-j\beta z}}={{E}_{0}}{{\operatorname{e}}^{-jkz}}+\frac{1}{2}jk\chi {{\operatorname{e}}^{-jkz}}\frac{C}{j(k-\beta )}-\frac{1}{2}jk\chi \frac{C}{j(k+\beta )}{{\operatorname{e}}^{-j\beta z}}-\frac{1}{2}jk\chi \frac{C}{j(k-\beta )}{{\operatorname{e}}^{-j\beta z}}$
We identify the term in ${{\operatorname{e}}^{-j\beta z}}$
$1=-\frac{1}{2}jk\chi \frac{1}{j(k+\beta )}-\frac{1}{2}jk\chi \frac{1}{j(k-\beta )}=-\chi \frac{{{k}^{2}}}{{{k}^{2}}-{{\beta }^{2}}}\to {{\beta }^{2}}={{n}^{2}}{{k}^{2}}={{n}^{2}}\frac{{{\omega }^{2}}}{{{c}^{2}}}$
So, we have $\beta =n\frac{\omega }{c}$
Identifying with the term in ${{\operatorname{e}}^{-jkz}}$ we find a second relation
$0={{E}_{0}}+\frac{1}{2}k\chi \frac{C}{(k-\beta )}$ . It gives $C=\frac{2(n-1)}{{{n}^{2}}-1}{{E}_{0}}=\frac{2}{n+1}{{E}_{0}}$
The field in the region $z>0$ is ${{E}_{t}}=\underbrace{\frac{2}{n+1}}_{t}{{E}_{0}}{{\operatorname{e}}^{j}}^{(\omega t-nkz)}$ and so we have the coefficient of transmission and the change of phase velocity by adding radiated fields with variable amplitudes, all propagating at velocity c.
But we can also compute the reflected waves in the $z<0$ region by adding all the radiated waves in this direction :
${{E}_{r}}(z)=-\frac{1}{2}jk\chi {{\operatorname{e}}^{jkz}}\int\limits_{0}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{-jkz'}}}=-{{\operatorname{e}}^{jkz}}\frac{{{n}^{2}}-1}{(n+1)}\frac{1}{n+1}{{E}_{0}}=-{{\operatorname{e}}^{jkz}}\frac{n-1}{n+1}{{E}_{0}}$
The reflected wave is ${{E}_{t}}(z,t)=\underbrace{\frac{1-n}{1+n}}_{r}{{E}_{0}}{{\operatorname{e}}^{j(\omega t+kz)}}$ and we have the coefficient of reflection without using any condition at the interface.
We can follow Feynman, chapter 31 of volume 1. The summary is :
A sinusoidal progressive plane wave in the vacuum ${{E}_{0}}{{e}^{i\omega (t-z/c)}}$ passes through an infinite plane of small thickness $\delta z$ and index $n$.
In the usual formalism of refraction index, it is delayed by $\delta t=(n-1)\delta z/c$ since it advances at the speed $c/n$ in the plate instead of $c$: time of travel in the plate $\delta z(n/c)$ instead of $\delta z(1/c)$.
The wave after the plate is $E={{E}_{0}}{{e}^{i(\omega (t-\delta t)-\omega z/c)}}={{e}^{i\omega (\delta t)}}{{E}_{0}}{{e}^{i(\omega (t-\delta t)-\omega z/c)}}=\underbrace{{{e}^{i\omega ((n-1)\delta z/c)}}}_{1-i\omega (n-1)\delta z/c}{{E}_{0}}{{e}^{i(\omega t-\omega z/c)}}$
Finally $E'={{E}_{0}}{{e}^{i(\omega t-z/c)}}\underbrace{-i\omega (n-1)\delta z/c{{E}_{0}}{{e}^{i(\omega t-z/c)}}}_{\text{Field radiated by the plane}}$
You just have to reverse the procedure: the sum of the incident wave and the radiated wave is indeed a delayed wave.
I recommend reading Feynman! (Sorry for my english)