function r=secant(f,x1,x2,geometry,dfdx,d2fdx2,varargin)
%SECANT Root finding, secant method.
%   X=SECANT(F,X1,X2) attempts to find the solution of F(X)=0 given 
%   X1 and X2 such that F(X1)*F(X2)<0 using the secant method. F is 
%   a function handle. 
%  
%   X=SECANT(F,X1,X2,GEOMETRY) the geometry of the method is ilustrated if
%   GEOMETRY='YES'.
%
%   X=SECANT(F,X1,X2,GEOMETRY,DFDX,D2FDX2) convegence plots are produced
%   where DFDX and D2FDX2 are handles for the first and second derivatives
%   of F.
%
%   X=SECANT(F,X1,X2,GEOMETRY,DFDX,D2FDX2,P1,P2,...) passes the parameters
%   P1,P2,... to F. If convergence analysis is not desired then set DFDX=[]
%   and D2FDX2=[].
%
%   See also NEWTON, BISECTION, IQI.

%   Francisco J. Beron-Vera, 2006/09/20
%   $Revision: 1.1 $  $Date: 2006/09/20 18:28:48 $

if nargin<4 
   geometry='no';
   ConvergenceAnalysis='no';
elseif nargin<5
   ConvergenceAnalysis='no';
else
   if isempty(dfdx) && isempty(d2fdx2)
      ConvergenceAnalysis='no';
   else
      ConvergenceAnalysis='yes';
   end
end

if f(x1,varargin{:})*f(x2,varargin{:})>=0
   error('F(X1)*F(X2) must be negative.')
end

if strcmp(geometry,'yes') || strcmp(geometry,'YES') % Plot function
   figure(1), clf
   X=linspace(x1,x2,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([x1 x2],[0 0],'k'), pause
   plot([x1 x1],[0 f(x1,varargin{:})],'g.:'), pause
   plot([x2 x2],[0 f(x2,varargin{:})],'g.:'), pause
   plot([x1 x2],f([x1 x2],varargin{:}),'r'), pause
end

% Secant loop.
N=24;
x=zeros(N,1);
x(1:2)=[x1 x2];
n=1;
while abs(x(n+1)-x(n))>eps*abs(x(n+1))
   n=n+1;
   if n>N
      warning('Tolerance not satisfied.')
      break
   end
   x(n+1)=...
      x(n)-...
      f(x(n),varargin{:})*(x(n)-x(n-1))/...
      (f(x(n),varargin{:})-f(x(n-1),varargin{:}));
   if strcmp(geometry,'yes') || strcmp(geometry,'YES')
      plot([x(n+1) x(n+1)],[0 f(x(n+1),varargin{:})],'g.:'), pause
      if f(x(n+1),varargin{:})*f(x(1),varargin{:})>0
         plot(x([2 n+1]),f(x([2 n+1]),varargin{:}),'r'), pause
      else
         plot(x([1 n+1]),f(x([1 n+1]),varargin{:}),'r'), pause
      end
   end

end
r=x(n);

if strcmp(ConvergenceAnalysis,'yes') % Convergence analysis.
   figure(2), clf
   x=x(1:n);
   n=0:n-1;
   e=abs(x-x(end));
   subplot(211)
   plot(n,e)
   xlabel('$n$',...
      'FontSize',14,...
      'Interpreter','LaTeX')
   ylabel('$e_n:=|x_n-r|$',...
      'FontSize',14,...
      'Interpreter','LaTeX')
   title(['$r=' sprintf('%16.8f',x(end)) '$'],...
      'FontSize',14,...
      'Interpreter','LaTeX')
   subplot(212)
   plot(n(2:end-1),e(2:end-1)./e(1:end-2).^1.62), hold on
   h=plot(n([2 end-1]),repmat(.5*d2fdx2(r)/dfdx(r),[1 2]),'r:'); hold off
   legend(h,{['$f''''(r)/2f''(r)\approx' num2str(.5*d2fdx2(r)/dfdx(r)) '$']},...
      'FontSize',14,...
      'Interpreter','LaTeX')
   legend boxoff
   xlabel('$k$',...
      'FontSize',14,...
      'Interpreter','LaTeX')
   ylabel('$e_n/e_{n-1}^{1.62}$',...
      'FontSize',14,...
      'Interpreter','LaTeX')
en