function w = wegstein(Func, x, tol, max)
%
%WEGSTEIN  Find a zero of a nonlinear equation of one variable 
%          using the wegstein algorithm. Func is the name of an 
%          M-file that contains a real-valued function of a  
%	   single real variable, and x is the initial guess.  
%          Note that a bad formulation of the problem may cause
%	   divergence.
%	   An optional third argument sets the tolerance for the 
%	   convergence. The presence of an optional fourth
%	   argument specifies the maximum iteration number.
%

%
% pre-processing, set default tolerance and maximum iteration  
% number.
%
if nargin<3, tol=0.000001; max=1000; end
if nargin==3, max=1000; end
t=[1:max+1]';

%
% initialization
%
x1 = x;
x2 = feval(Func, x1);

%
% start wegstein algorithm
%
for i=1:max,
    x3(i) =(x1*feval(Func,x2)-x2*feval(Func,x1))...
           /(x1-x2-feval(Func,x1)+feval(Func,x2));
    if(abs(x3(i)-x2) < tol)
	T = x3(i)
	iteration = i
 	t = [1:iteration+1]';
	break;
    end
    x1 = x2;
    x2 = x3(i);
end

%
% plot to see convergence/divergence
%
x3 = [x x3];
plot(t, x3, '-')
title('Convergence/Divergence of the Wegstein Algorithm')
xlabel('iteration number+1')

%
% examples of calling this function:
%
%	>> wegstein('cp2', 450)
%	>> wegstein('cp2', 450, 0.000001)
%	>> wegstein('cp2', 450, 0.000001, 1000)
%
% all these function calls in Matlab gives following results:
%
% 	T = 443.5712
%       iteration = 7
%
% together with a graph window.
