The Fresnel reflection for an arbitrary number of layers in the case of TM polarization is

\[\begin{align} \tilde{r}_{j,l,m \ldots n} = \frac{\tilde{r}_{j,l} + \tilde{r}_{l,m \ldots n} \exp\left(2 \mathrm{i} k_{z,l} d_l\right)} {1+\tilde{r}_{j,l} \tilde{r}_{l,m \ldots n} \exp\left(2 \mathrm{i} k_{z,l} d_l\right)} \end{align}\]

where \(d_l\) is the thickness of the \(l\)th metal layer and

\[\begin{align}\tilde{r}_{i,j}&= \left.\left(\frac{k_{z,i}}{\epsilon_i}-\frac{k_{z,j}}{\epsilon_j}\right) \middle/ \left(\frac{k_{z,i}}{\epsilon_i}+\frac{k_{z,j}}{\epsilon_j}\right)\right.\\ &=\frac{\epsilon_j k_{z,i}-\epsilon_i k_{z,j}} {\epsilon_j k_{z,i}+\epsilon_i k_{z,j}}\end{align}\] is the two layer Fresnel relation.

This equation can be iterated to produce the correct reflectivity for any number of layers, as the following Python snippet shows:

# n layer Fresnel
def r(kx,epsilon,d):
    # two layer Fresnel
    rtwo = lambda kx,epsilon: \
        (epsilon[1]*kz(kx,epsilon[0])-epsilon[0]*kz(kx,epsilon[1]))/ \

  if len(epsilon) == 2:
    return rtwo(kx,epsilon)

or in octave

% returns the complex n layer Fresnel reflectivity for TM polarization
function out = nlayerfresnel(k0,kx,epsilon,d)
    kz = @(kx,epsilon) sqrt(k0.^2.*epsilon-kx.^2);
    % two layer Fresnel
    r12 = @(kx,epsilon) \
                (epsilon(2).*kz(kx,epsilon(1))-epsilon(1).*kz(kx,epsilon(2)))./ \

    % if the length is two, return the two layer Fresnel, otherwise
    % recursively go through the layers
    if length(epsilon) == 2
        out = r12(kx,epsilon);
        out = (r12(kx,epsilon)+nlayerfresnel(k0,kx,epsilon(2:end),d(2:end)).*exp(2.0i.*kz(kx,epsilon(2)).*d(2)))./\

Note here that epsilon and d are ordered arrays holding the ordered list of complex permittivities of the system you’re trying to model.