function fern(N,MS,Kolor)

% This routine draws Barnsley's fern 
% The transformations for the fern are each in the form
% T(x,y) = (a*x+b*y+e, c*x+d*y+f). If the probability that a 
% particular transformation T is selected is p, you can characterize
% each transformation by providing the values of
% a,b,c,d,e,f, and p, in that order. 
%
% T1 -> 0,0,0,.16,0,0,.01
% T2-> .85,.04,-.04,.85,0,1.6,.85
% T3-> .2,-.26,.23,.22,0,1.6,.07
% T4-> -.15,.28,.26,.24,0,.44,.07
%
% INPUT
%	N -> The number of iterations.

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

% Initialize figure window and axes
fig1 = figure;
set(fig1,'Color','k','Position',[50   50   560   420])
ax1=axes;

% 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."
X(1,:)=[0.1,0.9];   

% The main loop where the iterations are performed.
for k=1:N-1
	r=rand;
	if r<.01
		X(k+1,:)=T(X(k,:),0,0,0,.16,0,0);
	elseif r<.86
		X(k+1,:)=T(X(k,:),.85,.04,-.04,.85,0,1.6);
	elseif r<.93
		X(k+1,:)=T(X(k,:),.2,-.26,.23,.22,0,1.6);
	else
		X(k+1,:)=T(X(k,:),-.15,.28,.26,.24,0,.44);
	end
end


plot(X(:,1),X(:,2),'o','markersize',MS, ...
    'MarkerFaceColor',Kolor,'MarkerEdgeColor',Kolor)

set(gca,'DataAspectRatio',[1 1 1])
% The transformation T

axis off

function U=T(X,a,b,c,d,e,f)
U=zeros(1,2);
U(1)=a*X(1)+b*X(2)+e;
U(2)=c*X(1)+d*X(2)+f;