Transfer matrix method for plasmonics

I run the shared the jreftran.m which has been designed for two layers:

delta=1i*g.*d.*ct;
M=zeros(2,2,length(d));
for j=1:length(d)
    M(1,1,j)=cos(delta(j));
    M(1,2,j)=1i./eta(j).*sin(delta(j));
    M(2,1,j)=1i*eta(j).*sin(delta(j));
    M(2,2,j)=cos(delta(j));
end    
M_t=[1,0;0,1]; %M total
for j=2:(length(d)-1)
    M_t=M_t*M(:,:,j);
end

But I intend to move to three or more layers using the same matrix method.