function [x1] = secant(f, a, b, tolerance, maxiter, printflag)
% secant  secant method
% f = function
% a = first initial guess
% b = second initial guess
% tolerance = maximum error desired
% maxiter = maximum iterations
% printflag = true, display intermediate results
% printflag = false, display only the final result

   x0 = a;
   fx0 = f(a);
   x1 = b;
   fx1 = f(b);
   
   % Iteration formula is
   % x1 <-- x1 - ( (x1 - x0)/(f(x1) - fx(0)) * f(x1)
   
   if printflag
      fprintf('%3s%24s%24s\n', 'n','x','err');
   end

   n = 1;
   while n <= maxiter
      dfx = fx1 - fx0;
      if dfx == 0
         error('division by zero');
      end
      
      dx = fx1*(x1 - x0) / dfx;
      err = abs(dx);
      
      x0 = x1;
      fx0 = fx1;
      x1 = x1 - dx;
      fx1 = f(x1);
      
      if printflag
         fprintf('%3d%24.15e%24.15e\n', n,x1,err);
      end
      
      if err < tolerance
         break;
      end
      
      n = n + 1;
   end
   if n > maxiter
      fprintf('Failure to converge in %d iterations\n', maxiter);
   end
end