function out=MultiNewton(F,x0,numits,tol)
%  Newton's Method defined for optimizing multi-dimensional F.
%
%     F- Name of a function that returns 
%        values of the function, the gradient AND the Hessian (in that
%        order)
%
%     x0= Initial value of x.
%
%     numits=  Maximum number of iterations
%
%     tol   =  Tolerance (like 1e-6) measures
%              change between iterates of x

for k=1:numits
    
   [g,gradg,hessg]=feval(F,x0);
   
   % This is an error check to see if the matrix is non-invertible
   if cond(hessg)>1000000
       error('The Hessian Matrix is not invertible');
   end
   
   xnew=x0-inv(hessg)*gradg;
   
   d=norm(xnew-x0);
   if d<tol
      out=xnew;
      break
   end
   
   x0=xnew;
end    
fprintf('Newton used %d iterations\n',k);
out=xnew;