S = [1/3 0; 0 1/3]; % the scaling matrix
d = sqrt(3)/3;
s = 2*sqrt(3)/9;
p_1 = [d; 0]; % each p_i is the coordinate of the bottom left corner of the six hexagons that remain after first step
p_2 = p_1 + 2*s*[1; 0]; % they are found using basic trigonometry; we assume the height of the outer pentagon is 1
p_3 = p_2 + 2*s*[1/2; sqrt(3)/2];
p_4 = p_3 + 2*s*[-1/2; sqrt(3)/2];
p_5 = p_4 + 2*s*[-1; 0];
p_6 = p_1 + 2*s*[-1/2; sqrt(3)/2];

F_1 = @(x) S*x + p_1; % the six functions to iterate
F_2 = @(x) S*x + p_2;
F_3 = @(x) S*x + p_3;
F_4 = @(x) S*x + p_4;
F_5 = @(x) S*x + p_5;
F_6 = @(x) S*x + p_6;

x0 = [1/2; 0]; % an arbitrary initial seed
n = 100000; % number of steps to iterate
orbit = zeros(2,n+1);
orbit(:,1) = x0;
for i = 1:n
    k = randi(6,1);
    if k == 1
        orbit(:,i+1) = F_1(orbit(:,i));
    elseif k == 2
        orbit(:,i+1) = F_2(orbit(:,i));
    elseif k == 3
        orbit(:,i+1) = F_3(orbit(:,i));
    elseif k == 4
        orbit(:,i+1) = F_4(orbit(:,i));
    elseif k == 5
        orbit(:,i+1) = F_5(orbit(:,i));
    else
        orbit(:,i+1) = F_6(orbit(:,i));
    end
end

scatter(orbit(1,:), orbit(2,:), 10, 'filled')