function [x,W,E]=TEA1(A,b,x0,n);
% TEA1 topological epsilon algorithm, first variant
%   [x,W,E]=TEA1(A,b,x0,n); topological epsilon algorithm to
%   accelerate the convergence of the sequence x_{j+1}=(I-A)*x_j+b
%   starting with x0. Returns the last result in x, the accelerated
%   vectors on the diagonal of the epsilon table in W, and the vectors
%   in the last row of the epsilon table in reverse order in E.

W=[]; r0=b-A*x0;
y=r0/norm(r0)^2; E(:,1)=x0;
for i=2:2*n+1
  v=zeros(size(x0));
  E(:,i)=(eye(size(A))-A)*E(:,i-1)+b;
  for j=i:-1:2
    if rem(i+j,2) == 0,                                     % i+j even
      de(:,j-1) = E(:,j)-E(:,j-1);
      w=v+y/(y'*de(:,j-1));
    else                                                    % i+j odd 
      w=v+de(:,j-1)/(de(:,j-1)'*(E(:,j)-E(:,j-1)));
    end;
    v=E(:,j-1);
    E(:,j-1)=w;
  end;
  W=[W w];
end
x=E(:,1);