function [yn,yns]=LambertNP(x)
% LAMBERTNP computes the inverse function yn 
%   of f(x)=x+exp(x) with Newton's method
%   and the derivative yns of the inverse function by
%   algorithmic differentiation
%   Uses the derivative fs(x)=1+exp(x) and fss(x)=exp(x)
tol=1e-6;
yas=0; ya=1; yns=0; yn=0;
while abs(ya-yn)>=tol
  yas=yns;ya=yn;
  yns=yas-(fs(ya)*(fs(ya)*yas-1)-fss(ya)*yas*(f(ya)-x))/fs(ya)^2;
  yn=ya-(f(ya)-x)/fs(ya);
end
yas=yns;ya=yn;         % one extra iteration 
yns=yas-(fs(ya)*(fs(ya)*yas-1)-fss(ya)*yas*(f(ya)-x))/fs(ya)^2;
yn=ya-(f(ya)-x)/fs(ya);