function [x, f, ok, count] = newtonND(fun, fungrad, funhess, xinit, epsi, maxiter)
if nargin < 6, maxiter = 30; end;
if nargin < 5, epsi = 100 * eps; end;
zero = 1e-8;
x0 = xinit;
count = 0;
while 1
   if count > maxiter
      break
   end;
   fhess = feval(funhess, x0);
   fgrad = feval(fungrad, x0);
   x = x0 - fhess \ fgrad;
   if max(abs(x - x0)) <= epsi * max(abs(x0))
      lambda = eig(fhess);
      if all(lambda > zero) 
         ok = 1;
      elseif all(lambda < -zero)
         ok = 2;
      else
         ok = 3;
      end;
      f = feval(fun, x);
      return
   end;
   x0 = x;
   count = count + 1;
end;
warning('przekroczona liczba iteracji w newtonND');
ok = 0;