function sierpinski(N)

% This routine draws a Sierpinski triangle with 
% vertices at A(0,0), B(1,0), and C(1,1).
% Note that the transformations T1, T2, T3 are written
% at the bottom of this routine as subfunctions. It is
% also assumed in this routine that each transformation
% has equal probability of being selected; i.e., p=1/3.
%
% INPUT
%	N -> The number of iterations.

% The vertices of triangle ABC.
A=[0,0];B=[1,0];C=[1,1];

% Close all existing figure windows. Delete this line
% if you do not want this action.
close all

% Initialize figure window and axes
fig1=figure;
ax1=axes;

% Plot and label vertices of triangle.
text(A(1),A(2),'A');
text(B(1),B(2),'B');
text(C(1),C(2),'C');

% Initialize some space to hold the sequence of points
% X1, X2, ..., XN
X=zeros(N,2);

% The initial value X1 = [0.5,0.5].
X(1,:)=[0.5,0.5];   % X(1,:) means "first row, every column."

% The main loop where the iterations are performed.
for k=1:N-1
	r=rand;
	if r<1/3
		X(k+1,:)=T1(X(k,:));
	elseif r<2/3
		X(k+1,:)=T2(X(k,:));
	else
		X(k+1,:)=T3(X(k,:));
	end
end

line(X(:,1),X(:,2),...
	'linestyle','none',...
   'marker','.',...
   'markersize',1)
	
% The transformations T1, T2, and T3.

function U=T1(X)
U=zeros(1,2);
U(1)=1/2*X(1);
U(2)=1/2*X(2);

function U=T2(X)
U=zeros(1,2);
U(1)=1/2*X(1)+1/2;
U(2)=1/2*X(2)+1/2;

function U=T3(X)
U=zeros(1,2);
U(1)=1/2*X(1)+1/2;
U(2)=1/2*X(2);