function y=LambertW(x)
% LAMBERTW implementation of the LambertW function
%   y=LambertW(x) computes the solution of x=y+exp(y) using
%   bisection

f=@(y) y+exp(y);
tol=1e-1;
a=-2; b=2;
y=(a+b)/2;
while (b-a)>tol,
  if f(y)<x, a=y; else b=y; end;
  y=(a+b)/2;
end;