function y=LambertN(x)
% LAMBERTN computes the inverse function 
%   of f(x)=x+exp(x) with Newton's method.
%   Uses the derivative fs(x)=1+exp(x)
tol=1e-12;
ya=1; yn=0;
while abs(ya-yn)>=tol
  ya=yn;
  yn=ya-(f(ya)-x)/fs(ya);
end
y=yn;