function [x,iter]=halley(p,x0,toll,nmax);
p1=polyder(p); p2=polyder(p1);
x(1)=x0; f=polyval(p,x(1)); fx=polyval(p1,x(1)); fxx=polyval(p2,x(1)); 
iter=1; test=1;
while test>toll &iter<=nmax
   iter=iter+1;
   rapp=f/fx;
   rapp1=(f*fxx)/(2*fx^2);
   f1=1/(1-rapp1);
   x(iter)=x(iter-1)-rapp*f1;
   test=abs(x(iter)-x(iter-1));
   f=polyval(p,x(iter)); fx=polyval(p1,x(iter)); fxx=polyval(p2,x(iter)); 
end