function [x,y,n]=mullerf(fn,xl,xr,xu,tol,max)
% Function to find the root of a function fn using the Muller method
%Inputs are:-
% fn = funtion name in quotes
% xl, xr, xu = 3 initial guesses of 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;
x(1)=xl;y(1)=feval(fn,xl);
x(2)=xr;y(2)=feval(fn,xr);
x(3)=xu;y(3)=feval(fn,xu);
i=1;
while error > tol
  c1=(y(i+1)-y(i))/(x(i+1)-x(i));
  c2=(y(i+2)-y(i+1))/(x(i+2)-x(i+1));
  d1=(c2-c1)/(x(i+2)-x(i));
  s=c2+d1*(x(i+2)-x(i+1));
  x(i+3)=x(i+2)-2*y(i+2)/(s+sign(s)*sqrt(s*s-4*d1*y(i+2)));
  y(i+3)=feval(fn,x(i+3));
  error=abs(y(i+3)); 
  n=n+1;i=i+1;
  if n >= max
      disp('Warning! Maximun number of iterations reached but not converged yet!')
      break
  end
end