You didn't explicitly set a value to R. If you had pre-allocated it, that would have prevented an R from a previous run to affect your data here. During debugging you can use clearvars, outside of debugging you should use functions.
R=zeros(size(theta));
for n=1:numel(theta)
lambda=632.8;%wavelength in vacuum in nm
k0=2*pi/lambda;
e1=1.30^2;e2=(0.15+1i*4)^2;e3=3.9^2;e4=1.33^2;e5=1.36;na=1.368;
n1=1.30;n2=sqrt(e2);n3=3.9;n4=1.33;n5=sqrt(e5);
d2=50;
d3=400;
d4=50;
theta=30:0.1:60;
x=n1*sind(theta);
q1=sqrt(e1-x.^2)./e1;
q2=sqrt(e2-x.^2)./e2;
q3=sqrt(e3-x.^2)./e3;
q4=sqrt(e4-x.^2)./e4;
b2=k0*d2.*q2;
b3=k0*d3.*q3;
b4=k0*d4.*q2;
m111=cosd(b2);
m112=-1i*sind(b2)./q2;
m121=-1i*q2.*sind(b2)./n1.^2;
m122=cosd(b2);
m211=cosd(b3);
m212=-1i*sind(b3)./q3;
m221=-1i*q3.*sind(b3);
m222=cosd(b3);
%//%M2=[m211,m212;m221,m222];
m311=cosd(b4);
m312=-1i*sind(b4)./q4;
m321=-1i*q4.*sind(b4);
m322=cosd(b4);
Mtot=zeros(1,length(theta));
for j=1:length(theta)
%M4=[m411(1,j),m412(1,j);m421(1,j),m422(1,j)];
M3=[m311(1,j),m312(1,j);m321(1,j),m322(1,j)];
M2=[m211(1,j),m212(1,j);m221(1,j),m222(1,j)];
M1=[m111(1,j),m112(1,j);m121(1,j),m122(1,j)];
Mtot=M3.*M2.*M1;
R(j) =((Mtot(1,1) +Mtot(1,2).*q4).*q1-Mtot(2,1) -Mtot(2,2).*q4)/((Mtot(1,1) +Mtot(1,2).*q4).*q1+Mtot(2,1) +Mtot(2,2).*q4).^2;
end
figure(1);
plot(theta,abs(R));
