function [TableOfSequence,TableOfIndexes]=...
    scaledvectorrho(S,indexes,method,toll,sequence)
% VECTOR RHO ALGORITHM FOR EXTRAPOLATION AT INFINITY OF SEQUENCE OF VECTORS
%
% --------------INPUT--------------------------------------------------------------
% S 		original sequence of N column vectors of length M (N=time) 
% indexes 	is a row vector of size N
% method	string= "originalsequence" or "integers" or "passedsequence"
% toll      avoid divisions 1/x with abs(x)<toll
% sequence	vector 1XN used only with method = "passedsequence" as  time series
%
% --------------OUTPUT-------------------------------------------------------------
% TableOfSequence	MXNXfloor(N/2)+1 table of values of accellerated sequence
% TableOfIndexes	corresponding time indexes

[M,N]=size(S);
% scaling of the sequence
Scaling=max(S(:,end),toll);
for j=1:N
    S(:,j)=S(:,j)./Scaling;
end

% 
TableOfSequence=zeros(N,M,N+2);TableOfSequence(:,:,2)=S';
TableOfIndexes=hankel(indexes,[indexes(end),zeros(1,length(indexes)-1)]);
switch method
    case 'originalsequence'
        for k=3:N+2
        for j=1:N-k+2%  -1?
           step=TableOfSequence(j+1,:,k-1)-TableOfSequence(j,:,k-1);
           norm2step=step*step';
           TableOfSequence(j,:,k)=...
               TableOfSequence(j+1,:,k-2) +(indexes(j+k-2)-indexes(j)).*...
               step./(norm2step);
        end
        end
    case 'integers'
        for k=3:N+2
        for j=1:N-k+2%  -1?
           step=TableOfSequence(j+1,:,k-1)-TableOfSequence(j,:,k-1);
           norm2step=step*step';
           TableOfSequence(j,:,k)=...
               TableOfSequence(j+1,:,k-2) +(k-2).*...
               step./(norm2step);
        end
        end
    case 'passedsequence'
        
        for k=3:N+2
        for j=1:N-k+2%  -1?
           step=TableOfSequence(j+1,:,k-1)-TableOfSequence(j,:,k-1);
           norm2step=step*step';
           TableOfSequence(j,:,k)=...
               TableOfSequence(j+1,:,k-2) +(sequence(j+k-2)-sequence(j)).*...
               step./(norm2step);
        end
        end
end
TableOfSequence=TableOfSequence(:,:,2:end);
% back-scaling of the table
for j=1:N
    for k=1:N
        vv=TableOfSequence(j,:,k);
        TableOfSequence(j,:,k)=vv.*Scaling';
    end;
end;
TableOfSequence=TableOfSequence(:,:,1:2:end);
TableOfIndexes=TableOfIndexes(:,1:2:end);

