%WOFWY The Lambert function w(y) for y >= 0
%
function w=wofy(y)
if y<0
  error('y out of range')
  return
end
if y<exp(1)
  w1=y;
else
  w1=log(y/log(y));
end
w0=-1;
it=0;
while abs(w1-w0)>10^2*eps
  it=it+1;
  w0=w1;
  w1=(w0^2+y*exp(-w0))/(1+w0);
end
w=w1;