function y=newton(f,dfdx,x,xmin,xmax,varargin)
%NEWTON Root finding, Newton method.
%   X=NEWTON(F,DFDX,X0,A,B) attempts to find the solution of F(X)=0 near 
%   X0 where DFDX is the derivative of the scalar function F defined 
%   in the interval [A,B] using Newton method. Both F and DFDX are function 
%   handles. The geometry of the method is ilustrated and convegence plots
%   are produced.
%
%   X=NEWTON(F,DFDX,X0,A,B,P1,P2,...) passes the parameters P1,P2,... to F.
%
%   See also BISECTION, SECANT, IQI.

%   Francisco J. Beron-Vera, 2006/09/07
%   $Revision: 1.1 $  $Date: 2006/09/07 09:00:25 $

% Plot function.
figure(1), clf
X=linspace(xmin,xmax,100);
plot(X,f(X,varargin{:}),'b')
axis tight, box off, hold on
xlabel('$x$',...
   'FontSize',14,...
   'Interpreter','LaTeX')
ylabel('$y$',...
   'FontSize',14,...
   'Interpreter','LaTeX')
disp('Press any key to continue.'), pause
plot([X(1) X(end)],[0 0],'k'), pause
plot([x x],[0 f(x,varargin{:})],'g.:'), pause

% Newton loop.
k=0;
y=x+1; Y=[]; 
while abs(y-x)>eps*abs(x)
   k=k+1;
   if k>50
      warning('Tolerance not satisfied.')
      break
   end
   y=x; Y=[Y y]; 
   x=x-f(x,varargin{:})/dfdx(x,varargin{:});
   plot(X,f(y,varargin{:})+dfdx(y,varargin{:})*(X-y),'r'), pause
   plot([x x],[0 f(x,varargin{:})],'g.:'), pause
end

% Convergence analysis.
figure(2), clf
k=0:length(Y)-1;
e=abs(Y-Y(end));
subplot(211)
plot(k,e)
xlabel('$n$',...
   'FontSize',14,...
   'Interpreter','LaTeX')
ylabel('$e_n:=|x_n-r|$',...
   'FontSize',14,...
   'Interpreter','LaTeX')
title(['$r=' sprintf('%16.8f',Y(end)) '$'],...
   'FontSize',14,...
   'Interpreter','LaTeX')
subplot(212)
plot(k(2:end-1),e(2:end-1)./e(1:end-2).^2)
xlabel('$k$',...
   'FontSize',14,...
   'Interpreter','LaTeX')
ylabel('$e_n/e_{n-1}^2$',...
   'FontSize',14,...
   'Interpreter','LaTeX')