function [x,y,n]=bisectf(fn,xl,xu,tol,max)
% Function to find the root of a function fn using the bisection method
%Inputs are:-
% fn = funtion name in quotes
% xl, xu =  lower and upper guesses which bracket the root
% tol = tolerance level
% max = maximum number of iterations allowed
%Outputs are:-
% x => desired root
% y => error in calculating root
% n => number of iterations
n=0;
error=1;
while error > tol
  fl=feval(fn,xl);
  fu=feval(fn,xu);
  xr=.5*(xu+xl);
  [xl,xu,xr];
  fr=feval(fn,xr);
  if fl*fr <= 0 
   xu=xr;
  elseif fr*fu <= 0
   xl=xr;
  else
   disp(' Error - root not bisected ')
   break
end
  x=xr;y=fr;
  error=abs(y); 
  n=n+1;
  if n >= max
      disp('Warning! Maximun number of iterations reached but not converged yet!')
      break
  end
end