function [y,it]=sroot(dig,p,pd,eps0,x0)
global s
digits(dig);
x=0; y=vpa(x0(s)); it=0;
while subs(abs(y-x))>eps0*subs(abs(y))
  it=it+1;
  if it>10
    s
    error('Newton diverges');
  end
  x=y;
  y=vpa(x-p(x,s)/pd(x,s));
end
