function [S1,n1,D,m1]=epsilonalgorithm(S,x,kmax)

% inizialization
R=cell(1,kmax);
n=cell(1,kmax);
R{1,1}=S;
N=size(S,1);
ns=repmat(x,1,size(R{1,1},2));
n{1,1}=ns;
% first new column
a=S(1:end-1);
d=diff(S);
d=diag(sum(d.^2,2).^(-1))*d;
%b=repmat(diff(x),1,size(R{1,1},2));
%R{1,2}=b.*d;
%R{1,2}=d.^(-1);
R{1,2}=d;
n{1,2}=ns(2:end,:);
% creating further columns
for k=2:kmax-1
    a=R{1,k-1};
    a=a(2:end-1,:);
    d=diff(R{1,k});
    d=diag(sum(d.^2,2).^(-1))*d;
    %R{1,k+1}=a+(d).^(-1);
    R{1,k+1}=a+d;
    newns=n{1,k};
    n{1,k+1}=newns(2:end,:);
end
j=0;
S1=cell(1,1);
n1=cell(1,1);
while 2*j+1<=size(R,2)
    S1{1,j+1}=R{1,2*j+1};
    n1{1,j+1}=n{1,2*j+1};
    j=j+1;
end
D=cell(1,length(x));
m1=cell(1,length(x));
h=1;
while h <= size(S1,2)
    mmax=size(S1{1,h},1);
    for m=1:mmax
        D{1,m}=[D{1,m};S1{1,h}(m,:)];
    end
    h=h+1;
end 
ncolumn=size(D{1,1},2);
 for i= 1:length(x)
     nrows=size(D{1,i},1);
     m1{1,i}=repmat((2.*(1:nrows)+i-1)',1,ncolumn);
     if mod(i,2)==0
         m1{1,i}=m1{1,i}-1;
     end
 end 
