Here is one way, where M_total_P is a 2x2x901 matrix, where 901 is the number of frequencies you have in the code.
clc
clear
close all
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n0 = 1; % Refractive index of air
ns = 1.46; % Refractive index of subtrate
% Parameters
lambda0=5e-6; % Wavelength of light in micrometers
%%% frequency Range
stsz = 0.01; %%% step size for frequency
frequency = 58:stsz:67; %% frequency in THz
c=3e8;n=25;
%%============= CHANGE HERE TO INITIALIZE M_total_P to size 2x2xlength(freq)
M_total_P = zeros(2,2,length(frequency));
%%=============
nA=1.685;dA = lambda0/(4*nA); %% Thickness of First Layer in meters
nB=1.501;dB = lambda0/(4*nB); %% Thickness of Second Layer in meters
for i = 1:length(frequency)
f = frequency(i);
w=2*pi*f*1e12; %%% Angular frequency by frequency
DAA=dA * nA * (w/c);DBB=dB * nB * (w/c);
%%% Transfer Matrix elements of first layer
ma11=cos(DAA); ma12=-1i*sin(DAA)/nA; ma21=-1i*nA*sin(DAA); ma22=cos(DAA);
MA=[ma11 ma12; ma21 ma22];
maa11=cos(DAA/2); maa12=-1i*sin(DAA/2)/nA; maa21=-1i*nA*sin(DAA/2); maa22=cos(DAA/2);
Ma=[maa11 maa12; maa21 maa22];
%%% Transfer MAtrix elements of Second layer
lb11=cos(DBB); lb12=-1i*sin(DBB)/nB; lb21=-1i*nB*sin(DBB); lb22=cos(DBB);
MB=[lb11 lb12; lb21 lb22];
lbb11=cos(DBB/2); lbb12=-1i*sin(DBB/2)/nB; lbb21=-1i*nB*sin(DBB/2); lbb22=cos(DBB/2);
Mb=[lbb11 lbb12; lbb21 lbb22];
%%============= CHANGE HERE TO STORE THE CORRECT VALUE IN CORRECT INDEX
M_total_P(:,:,i) = (Mb*MA*Mb)^n*(Ma*MB*Ma)^n;
%%=============
end
toc
size(M_total_P)
