function x = optn(objfun,gfun,x0,xtol,gtol)
% OPTN optimization by Newton's method with line search
% 
%   objfun  string, name of M-file that evaluates objective function
%   gfun    string, name of M-file that evaluates gradient+Hessian
%   x0      starting guess
%   xtol    tolerance for x
%   gtol    tolerance for gradient
%           Stop iterations when
%           || last correction in x || < xtol and
%           || gradient || < gtol
%   
npt = 25;            % # pts to show in linesearch plot
%  
x = x0;
f = feval(objfun,x);
[g,H] = feval(gfun,x);
dx = 1 + xtol;
while dx > xtol | norm(g) > gtol
% 
% Newton step
% 
  s = -(H\g);
  L = 1;
  xnew = x + L*s;
  fnew = feval(objfun,xnew);
  while fnew >= f
% 
%   Line Search
%   
    c = f; b = g'*s; a = (fnew-c-b*L)/L^2;
    lvals = linspace(0,L,npt);
    pvals = a*lvals.^2+b*lvals+c;
    fline = zeros(npt,1);
    for k=1:npt, fline(k)=feval(objfun,x+lvals(k)*s); end
    figure(2)
    plot(lvals,pvals,'r-', lvals,fline,'b*')
    xlabel('lambda'), ylabel('f(x+lambda*s)')
    title('Function along Newton direction (*) and parabolic approx')
    disp('Strike key to continue')
    pause
% 
    L = max(L/4,-b/(2*a));
    xnew = x + L*s;
    fnew = feval(objfun,xnew);
  end  
  [f fnew]
  dx = norm(L*s);
  figure(1)
  plot([x(1) xnew(1)],[x(2) xnew(2)],'r--')
  x = xnew;
  f = fnew;
  [g,H] =feval(gfun,x);
end