%SWOFY Symbolic version of wofy.m
%
function [w,it]=swofy(dig,y)
if subs(y)<0
  error('y out of range')
  return
end
eps0=.5*10^(-dig);
digits(dig+2);
if subs(y)<exp(1)
  w1=vpa(y);
else
  w1=log(vpa(y)/log(vpa(y)));
end
diff0=1;
it=0;
while abs(diff0)>eps0
  it=it+1;
  if it>200
    error('divergence in swofy')
    break
  end
  w0=w1;
  w1=vpa((w0^2+y*exp(-w0))/(1+w0));
  diff=abs((w1-w0)/w1);
  diff0=subs(diff);
end
w=w1;