function [yn,yns]=LambertW2P(x)
% LAMBERTW2P compute LambertW function and its derivative
%   y=LambertW2P(x) computes the solution of x=y+exp(y) using
%   Newton's method, and also the derivative of the LambertW
%   function using algorithmic differentiation

f=@(y) y+exp(y); fs=@(y) 1+exp(y); fss=@(y) exp(y);
tol=1e-6;                          % enough because of extra iteration
ya=1; yns=0; yn=0;                 % where Newton doubles accuracy    
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);
