 
w1=[0 0 0; 0 0.16 0; 0 0 1];
w2=[0.2 -0.26 0; 0.23 0.22 1.6; 0 0 1];
w3=[-0.15 0.28 0; 0.25 0.24 0.44; 0 0 1];
w4=[0.85 0.04 0; -0.04 0.85 1.6; 0 0 1];

% a=[0 0.2 -0.15 0.85];
% b=[0 -0.26 0.28 0.04];
% c=[0 0.23 0.25 -0.04];
% d=[0.16 0.22 0.24 0.85];

% e=[0 0 0 0]; % original code
% f=[0 1.6 0.44 1.6]; % original code
p=[0.01 0.09 0.16 1.0];

% step 1: translation
h=10; v=10;
L=[1 0 h; 0 1 v; 0 0 1];

% step 2: rotation
theta=pi/4;
R=[0.707 -0.707 0; 0.707 0.707 0; 0 0 1];

% step 3: dilation/contraction
bx=0.8; by=0.8;
D=[bx 0 0; 0 by 0; 0 0 1];

% step 4: shearing
kx=0.1; ky=0.1;
S=[1 kx 0; ky 1 0; 0 0 1];

% step 5: reflection

Fx=[-1 0 0; 0 1 0; 0 0 1]; % reflection with y axis.

CPS=L*R*D*S*Fx;

W1=CPS*w1*CPS^(-1);
W2=CPS*w2*CPS^(-1);
W3=CPS*w3*CPS^(-1);
W4=CPS*w4*CPS^(-1);


n=5000;            
            x=10; y=10;
                for l=1:1:n 
                pk=rand;
                 if pk <= p(1) 
                    K=W1;
                 elseif pk <= p(2)
                    K=W2;
                 elseif pk <= p(3)
                      K=W3;
                 else K=W4;
                 end
                 
                 xn(l) =  K(1,1)*x+ K(1,2)*y + K(1,3);
                 yn(l) =  K(2,1)*x+ K(2,2)*y + K(2,3);
                    x=xn(l);
                    y=yn(l);
                  
            end
            
            % subplot(2,2,2);
            plot(xn,yn,'.');
            % xlabel('(b)');
            axis equal;
            axis([-5 12 -1 20]);