Transfer matrix Method for trilayer

%TMM 3 Layers

% SOURCE PARAMETERS
c = 3e8;     % Speed of light in vacuum

% Operating frequencies
lmd = 800e-9;
freq = c./lmd;
omega = 2.*pi.*freq;

% Wavevector for free space
K0 = (2*pi)./lmd;

N = 1; % Number of layers excluding substrate and superstate

% Permivittity values
eps = [1,2.25,1];

% perambility values
mu = [1,1,1]; %permeability of every layer

% Refractive index
n = [sqrt((eps(1)*mu(1))),sqrt((eps(2)*mu(2))),sqrt((eps(3)*mu(3)))];

% Layer thickness
d = [600e-9,200e-9,600e-9];  %thickness of 2nd, 3rd and 4th layer in nm

% Angle of incidences
theta = 0:5:90;
thetaa = [0 0 0];   % Incident angles(in degrees)

% Wave vector
k2x = n(1,2).*(omega./c).*cosd(thetaa(1,2));

% Propagation Matrix
phi2 = k2x*d(1,2);
P2 = [exp(1i*phi2) 0; 0 exp(-1i*phi2)];

D1=zeros(2,length(theta));

% Dynamical matrix 1
D1_11 = 1;
D1_12 = 1;
D1_21 = n(1,1).*cosd(theta);
D1_22 = -n(1,1).*cosd(theta);
% D1 = [D1_11 D1_12; D1_21 D1_22]

% Dynamical matrix 2
D2_11 = 1;
D2_12 = 1;
D2_21 = n(1,2)*cosd(theta(1,2));
D2_22 = -n(1,2)*cosd(theta(1,2));
% D2 = [D2_11 D2_12; D2_21 D2_22]

% Dynamical matrix
D3_11 = 1;
D3_12 = 1;
D3_21 = n(1,3)*cosd(theta(1,3));
D3_22 = -n(1,3)*cosd(theta(1,3));
% D3 = [D3_11 D3_12; D3_21 D3_22]




for ii = 1:length(theta)
    D_1 = [D1_11(1,ii) D1_12(1,ii); D1_21(1,ii) D1_22(1,ii)]
    D_2 = [D2_11(1,ii) D2_12(1,ii); D2_21(1,ii) D2_22(1,ii)]
    D_3 = [D3_11(1,ii) D3_12(1,ii); D3_21(1,ii) D3_22(1,ii)]
    
    M = D_1^(-1)*D_2*P2*D_2^(-1)*D_3;
    
    rs = M(2,1)./M(1,1)
    ts = 1./M(1,1)
    Rs = abs(rs.^2)
    Ts = ((n(1,3).*cosd(thetaa(1,3))./(n(1,1).*cosd(theta))).*abs(ts.^2))
    Frensel_s = Rs + Ts
end

D_1 = 2×2

1 1 1 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D_2 = 2×2

1.0000 1.0000 1.4943 -1.4943

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D_3 = 2×2

1.0000 1.0000 0.9848 -0.9848

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rs = -0.2062 + 0.1971i

ts = -0.6546 - 0.7102i

Rs = 0.0814

Ts = 1×19

0.9328 0.9364 0.9472 0.9657 0.9927 1.0292 1.0771 1.1387 1.2177 1.3192 1.4512 1.6263 1.8656 2.2072 2.7273 3.6041 5.3718 10.7027 Inf

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Frensel_s = 1×19

1.0142 1.0177 1.0286 1.0471 1.0740 1.1106 1.1585 1.2201 1.2991 1.4005 1.5326 1.7077 1.9470 2.2886 2.8087 3.6854 5.4532 10.7841 Inf

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Index in position 2 exceeds array bounds. Index must not exceed 1.

% Dynamical Matrix for S - Wave

plot(theta,Rs,theta,Ts)
legend('R','T')

%% Graph output should be like this

[Walter Roberson]

D1_11 = 1;

That is a 1 x 1 scalar.

for ii = 1:length(theta)
    D_1 = [D1_11(1,ii) D1_12(1,ii); D1_21(1,ii) D1_22(1,ii)]

That attempts to index the 1 x 1 scalar at position (1,ii) which is going to fail when ii is 2.