aa=0.02;  % set angle increment. You may pick larger aa to 
          % get fewer points 

th=-pi:aa:pi;
N=length(th);
x1=zeros(2,N);
x2=x1;
x8=x1;
x4=x1;
for i=1:N;
    a=sin(th(i));b=cos(th(i)); % x0=[a;b];
    c1=abs(a)+abs(b);  % one norm of x0
    c8=max([abs(a) abs(b)]); % infinity norm of x0
    c3=(abs(a)^3+abs(b)^3)^(1/3); % three norm
    x1(:,i)=[a;b]/c1;  % set the i-th column of x1 
    x8(:,i)=[a;b]/c8;
    x3(:,i)=[a;b]/c3;
    x2(:,i)=[a;b];  % no need to scale since the two norm of [a;b] is 1.
end

S1=[1 -1;1 0];  % Your S1 is different
s1x1=S1*x1;s1x8=S1*x8; s1x3=S1*x3; s1x2=S1*x2;
%  You need to do the same for S2 and S3;

% The following only plots two rows of figures
% You need to use two more rows for the maps under S2 and S3
%  and use subplot(4,4,n), n from 1 to 16

subplot(2,4,1)
plot(x1(1,:),x1(2,:),x1(1,:)*2,x1(2,:)*2,x1(1,:)*3,x1(2,:)*3)

subplot(2,4,2)
plot(x2(1,:),x2(2,:),x2(1,:)*2,x2(2,:)*2,x2(1,:)*3,x2(2,:)*3)

subplot(2,4,3)
plot(x3(1,:),x3(2,:),x3(1,:)*2,x3(2,:)*2,x3(1,:)*3,x3(2,:)*3)

subplot(2,4,4)
plot(x8(1,:),x8(2,:),x8(1,:)*2,x8(2,:)*2,x8(1,:)*3,x8(2,:)*3)

subplot(2,4,5)
plot(s1x1(1,:),s1x1(2,:),s1x1(1,:)*2,s1x1(2,:)*2,s1x1(1,:)*3,s1x1(2,:)*3)

subplot(2,4,6)
plot(s1x2(1,:),s1x2(2,:),s1x2(1,:)*2,s1x2(2,:)*2,s1x2(1,:)*3,s1x2(2,:)*3)

subplot(2,4,7)
plot(s1x3(1,:),s1x3(2,:),s1x3(1,:)*2,s1x3(2,:)*2,s1x3(1,:)*3,s1x3(2,:)*3)

subplot(2,4,8)
plot(s1x8(1,:),s1x8(2,:),s1x8(1,:)*2,s1x8(2,:)*2,s1x8(1,:)*3,s1x8(2,:)*3)