l0=729;
a=[];a(1)=0;a(2)=l0;
x0=a(1):a(2);
y0=0.*ones(1,730);
plot(x0,y0,'r');
hold on;
for j=1:5
l=l0/3^j;
p1=1:4:2^(j+1);
p2=2:4:2^(j+1);
p3=3:4:2^(j+1);
p4=4:4:2^(j+2);
A=[];
for i=1:2^(j+1)
p11=fix((i+3)/4);p22=fix((i+2)/4);p33=fix((i+1)/4);p44=fix(i/4);
if i==p1(p11)
A(i)=a((i+1)/2);
elseif i==p2(p22)
A(i)=a(i/2)+l;
elseif i==p3(p33)
A(i)=a((i+1)/2)-l;
elseif i==p4(p44)
A(i)=a(i/2);
end;
end;
for m=1:2^(j+1)
if mod(m,2)==1
xm=A(m):A(m+1);
ym=j.*ones(1,A(m+1)-A(m)+1);
plot(xm,ym,'r');
hold on;
a(m)=A(m);
a(m+1)=A(m+1);
end;
end;
end;
xlim([-20,750]);
ylim([-1,j+1])