function Lnf = newton(x0,f0,x)
  d = divdiff(x0,f0);
  ds = d(1,:);
  n = length(ds);
  m = length(x);
  Lnf = zeros(size(x));
  for i=1:m
    for j=1:n
      Lnf(i) =  Lnf(i) + ds(j) * (prod( x(i) - x0([1:j-1])));
    endfor

  endfor
end;