clear

p=[0 0;10 0];n=1;

A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)];

for k=1:5

   j=0;

   for i=1:n

      q1=p(i,:);

      q2=p(i+1,:);

      d=(q2-q1)/3;

      j=j+1;r(j,:)=q1;

      j=j+1;r(j,:)=q1+d;

      j=j+1;r(j,:)=q1+d+d*A';

      j=j+1;r(j,:)=q1+2*d;

   end

   n=4*n;clear p

   p=[r;q2];

end

plot(p(:,1),p(:,2))