% fractal through iterated random functions

N=10000;
x=zeros(2,N);
w=zeros(2,1);
w=[0.2993; 0.7007];
A=zeros(2,2,2);
b=zeros(2,2);

A(:,:,1)=[0.4 -0.3733; 0.06 0.6];
b(:,1)=[0.3533; 0.0];
A(:,:,2)=[-0.8 -0.1867; 0.1371 0.8];
b(:,2)=[1.1; 0.1];

u=rand(N,1);

x(:,1)=rand(2,1)-0.5*ones(2,1);

for i=2:N
    if (u(i,1)<w(1,1))
        x(:,i)=A(:,:,1)*x(:,i-1)+b(:,1);
    else
        x(:,i)=A(:,:,2)*x(:,i-1)+b(:,2);
    end
end
plot(x(1,:),x(2,:),'.')