function [Y,X,U,Gamma]=MPE(G,d,x0,n)
% MPE Minimal Polynomial Extrapolation 
%   [Y,X,U,Gamma]=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, and Gamma contains the coefficients of the
%   polynomial at each step.
     
n=length(d);
Y=[]; U=[]; Gamma=zeros(n+1,n+1);
xk=x0; X=xk;
for k=0:n
  xk1=G*xk+d; X=[X,xk1];
  uk=xk1-xk; U=[U,uk];
  c=[-U(:,1:k)\uk;1]; 
  y=X(:,1:k+1)*c/sum(c);
  Gamma(:,k+1)=[c/sum(c); zeros(n-k,1)];
  Y=[Y,y];
  xk=xk1;
end;  