function y=Lambert(x)
% LAMBERT computes the inverse function of 
%    f(x)=x+exp(x) by bisection
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;