function [Y,X,U]=TEA(G,d,x0,q,n)
% TEA Topological Epsilon Algorithm  
%   [Y,X,U]=mpe(G,d,x0,n) performs n number of iteration
%   steps x(j+1)=G * x(j) + d starting with x0 and computes an
%   extrapolated vector y at each iteration, stored in the matrix
%   Y. The basic iteration steps and their differences are stored
%   in X and U respectively.

X=[x0]; Y=[x0]; U=[]; xk=x0; QU=1;
for k=1:n
  m=k+1;                   % size of linear system 
  xk1=G*xk+d; X=[X,xk1];   % compute two basic iterations  
  xk2=G*xk1+d; X=[X,xk2];
  uk1=xk1-xk; uk2=xk2-xk1; % form the new differences
  U=[U,uk1, uk2];
  QU(1,m)=1;               % update the matrix for gamma
  for j=2:m-2              % copy elements in new column m
    QU(j,m)=QU(j+1,m-1);
  end
  for j=1:m-2              % copy elements for new row m
    QU(m,j)=QU(m-1,j+1);
  end
  QU(m,m-1)=q'*uk1;        % update elements in corner 
  if m>2, 
    QU(m-1,m)=QU(m,m-1); 
  end
  QU(m,m)=q'*uk2;
  ee=zeros(m,1); ee(1)=1;  % right hand side 
  gamma=QU\ee;
  y=X(:,1:k+1)*gamma;      
  Y=[Y,y];
  xk=xk2;
end;  