N=10000;
p=4;
a=zeros(2,1,4);
a(:,:,1)=[1 ; 1];
a(:,:,2)=[-1 ; 1];
a(:,:,3)=[-1 ; -1];
a(:,:,4)=[1 ; -1];

A=diag([1/3,1/3]);
x=zeros(2,1);

for i=1:N
    x=A*x+ a(:,:, floor(1+p*rand(1,1,'def')))*2/3  ;
    plot2d(x(1),x(2),0,rect=[-1.1,-1.1,1.1,1.1]);
end