clc;
clear;

x = [0:0.1:1]

for j = 1:length(x)

         y(j) = (j-1)*10;

        m11(j) = cos(y(j));
        m12(j) = -sin(y(j));
        m21(j) = sin(y(j));
        m22(j) = cos(y(j));

       M = [m11(j) m12(j);m21(j) m22(j)];

       n11(j) = -cos(y(j));
       n12(j) = sin(y(j));
       n21(j) = sin(y(j));
       n22(j) = -cos(y(j));

       N = [n11(j) n12(j);n21(j) n22(j)];

       MN= M.*N;

       P= MN.^8

       p11(j) = P(1,1);
       p12(j) = P(1,2);
       p21(j) = P(2,1);
       p22(j) = P(2,2);

      R(j) = (p11(j)+p22(j))/(p21(j)+p22(j));
end

plot(y,R)