function val = WynnEpsilon(a,steps,object,display)

% function S = WynnEpsilon(a,steps,object,display) 
%
% extrapolates a sequence or series by using Wynn's epsilon algorithm
%
%   a           vector of terms
%   steps       number of extrapolation steps
%   object      'sequence': extrapolation to the limit of a
%               'series'  : extrapolation to the limit of cumsum(a)
%   display     'on':   plot the absolute value of differences in the
%               first row of the extrapolation table for diagnosis
%
%   val           extrapolated value

S = zeros(length(a),2*steps+1);
switch object
    case 'sequence'
        S(:,1) = a;
    case 'series'
        S(:,1) = cumsum(a);
    otherwise
        error('MATLAB:badopt','%s: no such object known',object);
end
for j=2:2*steps+1
    if j==2
        switch object
            case 'sequence'
                S(1:end-j+1,j) = 1./diff(a);
            case 'series'
                S(1:end-j+1,j) = 1./a(2:end);
        end
    else
        S(1:end-j+1,j) = S(2:end-j+2,j-2) + 1./diff(S(1:end-j+2,j-1));
    end
end
S = S(:,1:2:end);
if isequal(display,'on')
    h=semilogy(0:length(S(1,:))-2,abs(diff(S(1,:))));
    set(h,'LineWidth',3);
end
val = S(1,end);

return


