function [x,H]=newton(fun,x0,opts,varargin)
  terminate=0;
  numiter=0;
  x=x0;
%  fprintf('starting point: '); x
  
  eps_fun=optimget(opts,'TolFun',1e-8);
  disp=optimget(opts,'Display','off');
  maxnumiter=optimget(opts,'MaxIter',100);
  bt_alpha=0.2; % (0,0.5) smaller alpha meaning easier to take a step
  bt_beta=0.8;  % (0,1) smaller beta meaning taking a smaller step
  while ~terminate,
    numiter=numiter+1;

    [f,g,H]=feval(fun,x,varargin{:});
    dx=-H\g;
    lambda_sq=(-g'*dx);
    if (numiter>maxnumiter)
      terminate=1;
    elseif (lambda_sq/2<eps_fun) %% lambda_sq/2=g'*H^-1*g/2 .= |min(parabola)-f(x)|
      terminate=1;
    else
      %% backtracking line search
      t=1;
      f_t=feval(fun,x+t*dx,varargin{:});
      while f_t>f+bt_alpha*t*(-lambda_sq),
        t=bt_beta*t;
        f_t=feval(fun,x+t*dx,varargin{:});
      end
      x=x+t*dx;
    end

    if strcmp(disp,'iter')
      fprintf('iter: %i, norm of g: %.1e, hesnorm of g: %.1e, stepsize: %g\n',numiter,norm(g),lambda_sq,t);
    end
  end
  
  % returned H