function [x,dx,err,it] = NEWTON(gradf,Hessef,x0,dxmax,gradfmax,itmax)

% [x,dx,err,it] = NEWTON(gradf,Hessef,x0,dxmax,gradfmax,itmax)
% 
% Newton-Verfahren zur Minimierung einer C^2-Funktion f
%
% Input
% -----
% gradf         Name der Funktion, die den Gradienten liefert
%
% Hessef        Name der Funktion, die die Hessematrix liefert
%
% x0            Startpunkt (als Spaltenvektor)
%
% dxmax         Abbruchkriterium Schrittweite
% gradfmax      Abbruchkriterium Norm Gradient
% itmax         Abbruchkriterium Anzahl Iterationen
%
% iteriere, solange norm(x(:,k+1) - x(:,k)) > dxmax
%               und norm(gradf(x(:,k+1)) > gradfmax
%               und k < itmax
%
% Output
% ------
% x             Matrix der Iterationspunkte 
%               x(:,k) ist Punkt vor k-ter Iteration
%
% dx            Vektor der Schrittweiten
%
% err           Fehlervektor, err(k) = norm(gradf(x(:,k))
%
% it            1 + Anzahl ausgeführter Iterationen


% Inputstrings --> Funktionen
gradf = str2func(gradf);
Hessef = str2func(Hessef);

% Startwerte
x=x0;
gf = gradf(x(:,1));
Hf = Hessef(x(:,1));
dx(1) = inf;            % dx > dxmax in jedem Falle erfüllt
err(1) = norm(gf);     
it = 1;

while dx(it) > dxmax && err(it) > gradfmax && it <= itmax
    
    % neue Iteration
    it = it+1;    
    
    % Löse NEWTON-Gleichung    
    opts.SYM = true;                            % Hessematrix ist symmetrisch
    [p, condition] = linsolve(Hf,-gf,opts);     % berechne Schrittrichtung
    
    % Abbruch, wenn Hessematrix schlecht konditioniert oder nicht invertierbar
    if condition < eps
        err(it) = inf;
        return
    end % if
    
    % neuer Iterationspunkt
    x(:,it) = x(:,it-1)+p;    
    
    % neue Daten
    gf = gradf(x(:,it));
    Hf = Hessef(x(:,it));
    
    % Abbruchkriterien
    dx(it) = norm(x(:,it)-x(:,it-1));
    err(it) = norm(gf);
    
end % while

it = it-1;