%This function calculates an affine Reccurent Fractal Interpolation Function.
%Input:IP, AD, S, N, steps
%IP is a (N+1)x2 matrix representing the interpolation points
%N is the number of mappings and N+1 is the number of points
%S is a Nx1 matrix representing the contraction factors
%AD: is a Nx4 matrix representing the interval that is mapped to each section.
%i.e. if [...;x_0,x_1,y_0,y_1;...] then the [x_0,x_1] is the interval.
%steps: is the number of iterations
%Output: P
%The algorithm computes and returns N^(steps+1)+1 points
%representing the attractor of the RFIF.
function P=myRFIF(IP,AD,S,N,steps)
%Here are the map parameters.
for i=1:N
a(i)=(IP(i+1,1)-IP(i,1))/(AD(i,2)-AD(i,1));
e(i)=(AD(i,2)*IP(i,1)-AD(i,1)*IP(i+1,1))/(AD(i,2)-AD(i,1));
c(i)=(IP(i+1,2)-IP(i,2))/(AD(i,2)-AD(i,1))-S(i)*(AD(i,4)-AD(i,3))/(AD(i,2)-AD(i,1));
f(i)=(AD(i,2)*IP(i,2)-AD(i,1)*IP(i+1,2))/(AD(i,2)-AD(i,1))-S(i)*(AD(i,2)*AD(i,3)-AD(i,1)*AD(i,4))/(AD(i,2)-AD(i,1));
end;
A_0=IP;
points=N+1;
E=10^-14;
for i=1:steps
%for each section
l=1;
for j=1:N
%find the start of the respective interval
for start=1:points
if abs(A_0(start,1)-AD(j,1))