% [xmin,fxmin]=plainnewton(f,gradf,hessf,x0,tol)
%
% This is plain old Newton's method used to minimize func.  
% The basic algorithm is   x^{i+1} = x^i - inv(hessf(x_i))gradf(x^i)
%
function [xmin,fxmin]=plainnewton(func,grad,hess,x0,tol)


%some useful constants
  iterlim = 200;   % maximum number of iterations allowed

 
%
%  Initialize x and oldx.  We use oldx=1000*x just to make sure that
%  it isn't too close to x.  
%
x=x0;
fx=func(x);
oldx=1000*x0;
oldfx=fx*1000;
%
% The main loop.  While the current solution isn't good enough, keep
% trying...  Stop after 200 iterations in the worst case.
%
  iter=0;
  tic; 
  df = grad(x);  % Compute the gradient
  ndf = sqrt(df'*df);   % Compute the norm of the gradient

 % Set things up to display each iteration
  disp('   iter         f(x)              ||df||');
  s=sprintf('     %d        %7.8f       %7.8f           %7.5f', iter, fx, ndf);
 disp(s)
  while (((abs(fx-oldfx) > tol*(1+abs(fx))) & (ndf > tol^2*(1+abs(fx)))) & (iter < iterlim))
%
%  Compute the hessian.
%
    H=hess(x);
%
%  Compute the rhs=-grad f(x)
%
    rhs=-df;
%
%  *****************************************************************************
%  *****************************************************************************
%  solve the equations to find  search direction d one of two ways, note you 
%  should have exactly one line commented out to assign d
% 
%  one way to do this is to solve  d = inv(H)*rhs which involves calculating 
%  the inverse of the Hessian.
     d = inv(H)*rhs;
%
%  another way to do this is to use Guassian elimination to solve the system 
%  of linear equations.
%      d=H\rhs;  
%  *****************************************************************************
%  *****************************************************************************
%
%   compute new point
    xnew = x +d;
%  The new point is at xnew.  Update oldx and other  info.
%

    oldx=x;
    oldfx=fx;
    x=xnew;
    fnew = func(x);
    fx = fnew;
    df = grad(x);

     
%  Display intermediate results.
%
  ndf = sqrt(df'*df);   % Compute the norm of the gradient

  iter=iter+1;
  % Display each iteration
  s=sprintf('     %d        %7.8f       %7.8f         %7.5f', iter, fx, ndf);
  disp(s);
%
%  end of the loop.
%
  end;
%
%  Setup the return values.
%
xmin=x;
fxmin=fx;

%
%  Print out some information
%
toc
disp('Done with Newtons method');