function [xo,Ot,nS]=newtont(S,x0,ip,G,H,Lb,Ub,problem,tol,mxit)
//   Unconstrained optimization using Newton with trust region.
//
//   [xo,Ot,nS]=newtont(S,x0,ip,G,H,Lb,Ub,problem,tol,mxit)
//
//   S: objective function
//   x0: initial point
//   ip: (0) no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
//   G: gradient vector function
//   H: Hessian matrix function (default = sparse finite-differences)
//   Lb, Ub: lower and upper bound vectors to plot (default = x0*(1+/-2))
//   problem: (-1): minimum (default), (1): maximum
//   tol: tolerance (default = 1e-4)
//   mxit: maximum number of iterations (default = 50*(1+4*~(ip>0)))
//   xo: optimal point
//   Ot: optimal value of S
//   nS: number of objective function evaluations

//   Copyright (c) 2001 by LASIM-DEQUI-UFRGS
//   $Revision: 1.0 $  $Date: 2001/07/16 21:14:02 $
//   Argimiro R. Secchi (arge@enq.ufrgs.br)

 if nargin < 4,
   error('newtont requires four input arguments!');
 end
  if isempty(ip),
   ip=0;
  end
 if isempty(G),
   error('newtont: gradient function must be provided!');
 end
 if nargin < 5 || isempty(H),
   hes=0;
 else
   hes=1;
 end
 if nargin < 6 || isempty(Lb),
   Lb=-x0-~x0;
 end
 if nargin < 7 || isempty(Ub),
   Ub=2*x0+~x0;
 end
 if nargin < 8 || isempty(problem),
   problem=-1;
 end
 if nargin < 9 || isempty(tol),
   tol=1e-4;
 end
 if nargin < 10 || isempty(mxit),
   mxit=50*(1+4*~(ip>0));
 end

 x0=x0(:);
 n=size(x0,1); 

 if ip && n == 2,
   global X_PLOT
   figure(abs(ip));
   [X1,X2]=meshgrid(Lb(1):(Ub(1)-Lb(1))/20:Ub(1),Lb(2):(Ub(2)-Lb(2))/20:Ub(2));
   [n1,n2]=size(X1);
   f=zeros(n1,n2);
   for i=1:n1,
    for j=1:n2,
      f(i,j)=evstr(S+"([X1(i,j);X2(i,j)])");
    end
   end
   mxf=max(max(f));
   mnf=min(min(f));
   df=mnf+(mxf-mnf)*(2.^(([0:10]/10).^2)-1);
   contour(X1,X2,f,df); set(gca(),"auto_clear","off");
//   clabel(v,h);
//   h1=plot(x0(1),x0(2),'ro');
   xplot(x0,0); h1=gca();
   legend(h1,'start point');
   X_PLOT=1;
   
   if ip > 0,
     disp('Pause: hit any key to continue...'); halt;
   end
 end
 
 if problem == 1,
   pr="-evstr(";
 else
   pr="evstr(";
 end;
 
   op=optimset('LargeScale','on','GradObj','on','MaxFunEvals',inf,'MaxIter',mxit,...
               'Display','off','TolFun',tol,'TolX',tol);
  
 if ip && n == 2,
     fun={pr+S+"(x);xplot(x,"+int2str(ip)+")*[]",pr+G+"(x)"};
 else
     fun={pr+S+"(x)",pr+G+"(x)"};
 end

 if hes,
     op=optimset(op,'Hessian','on');
     fun=[fun,{pr+H+"(x)"}];
 end

 [xo,yo,exitflag,out] = fminunc1(fun,x0,op);
 nS=out.funcCount;
  it=out.iterations;

 Ot=-yo*problem;
 
 if it >= mxit,
   disp(['Warning Newtont: reached maximum number of iterations!']);
 elseif ip && n == 2,
   figure(abs(ip)); plot(xo(1),xo(2),'r*'); h2=gca();
   legend([h1,h2],'start point','optimum');
 end
