

function levy(n)

%levy(16),iteration=16
%initial conditions
n=1;

x1=0;y1=0;

x2=1;y2=0;

%set up vetors

[X,Y]=levy1(x1,y1,x2,y2);


for i=1:n

    for j=1:length(X)/3
        
        w=levy2(X(1+3*(j-1):3*j),Y(1+3*(j-1):3*j));

        [XX(3*2*(j-1)+1:3*2*(j-1)+3),YY(3*2*(j-1)+1:3*2*(j-1)+3)]=levy1(w(1,1),w(1,2),w(1,3),w(1,4));

        [XX(3*2*(j-1)+3+1:3*2*(j-1)+3+3),YY(3*2*(j-1)+3+1:3*2*(j-1)+3+3)]=levy1(w(2,1),w(2,2),w(2,3),w(2,4));
    end

    X=XX;
    Y=YY;
end

plot(X,Y)

hold on

axis equal

 
% find the midpoint of (x1,y1),(x2,y2),then save the point into the vector

function [X,Y]=levy1(x1,y1,x2,y2)

x3=1/2*(x1+x2+y1-y2);

y3=1/2*(-x1+x2+y1+y2);

X=[x1,x3,x2];

Y=[y1,y3,y2];

 

%divide x,y points into 2 groups
% each represents the end point of a line
% save them it w, which represents the location of of iteration n. 

function w=levy2(X,Y)

a11=X(1);b11=Y(1);

a12=X(2);b12=Y(2);

a21=X(2);b21=Y(2);

a22=X(3);b22=Y(3);

w=[a11,b11,a12,b12;a21,b21,a22,b22];

