x=[0;0];
N=5000;

figure(1) 
axis([-5 5 0 12]);

plot(x(1),x(2),'.')

hold on
% Enter each affine transformation Ax+T
A1=[0 0; 0 .16];
T1=[0;0];
A2=[.85 .04; -.04 .85];
T2=[0;1.6];
A3=[0.2 -.26; .23 .22];
T3=[0;1.6];
A4=[-.15 .28; .26 .24];
T4=[0;.44];
% Enter the probabilities
Pr=[.01,.85,.07,.07];
% The iteration
for i=1:N
    r = rand;
    if r < Pr(1) 
       y = A1*x+T1;    
    elseif  r < Pr(1) + Pr(2) 
       y = A2*x+T2; 
    elseif r < Pr(1) + Pr(2)+Pr(3) 
       y = A3*x+T3; 
    else
       y = A4*x+T4; 
    end 
    x=y;
    plot(x(1),x(2),'.')
    drawnow
 end
hold off