function [x,fehler,k] = my_Newton(f,x,h,n,err)

% Die Funktion my_Newton.m liest eine andere Funktion f: R^n -> R^2 ein, 
% berechnet die Jacobi Matrix von f und fuehrt das Newton-Verfahren zur
% Bestimmung von Nullstellen durch. Dabei wird das LGS:

%         f'(x) d = -f'(x)

% mit dem Gaussverfahren geloest. 

% benoetigte Variablen:
% f         Zeiger auf diese Funktion f
% x         Startwert der Newton-Iterationen
% h         Schrittweite
% n         Raumdimension
% err       geforderte Genauigkeit des Newton-Verfahrens

% Ausgaben:
% x         Nullstelle
% fehler    Fehler ||f(x)|| in der 2-Norm
% k         Anzahl der benoetigten Iterationen 



% Bilde Einheitsmatrix fuer die Berechnung des zentralen Differenzen-
% quotienten
E  = eye(n);

% Startwert, 1. Iteration
k=1;

% Fehler in den Startwerten:
fehler(1)=norm(f(x),2);


while fehler(k) > err
    
    % rechte Seite des LGS:
    b = - f(x);   
    
    % Bilde Jacobi-Matrix (von Hand): 
    A=zeros(n,n);
    for i=1:n
        for j=1:n
            l = f(x + h*E(:,j)); % E(:,j): die j.Spalte von E
            r = f(x - h*E(:,j)); % l und r sind Spaltenvektoren der Laenge n
            A(i,j)= l(i)-r(i);
        end
    end
    A = A/(2*h);

    % Verwende die effiziente Matlab Routine "\" zur Loesung des LGS:
    % d = A\b   
    
    % Verwende die selbst programmierte Gausselimination aus Aufgabe 8: 
    % (mit Pivotsuche, falls pivot = 1, anderenfalls ohne)
    index = [1:n]';
    pivot = 1; 
    [A,d,index] = gaussElim(A,b,n,index,pivot); 
    
    % Ueberschreibe das alte x mit dem neuen Iterationswert: 
    x = x + d;

    % Erhoehe die Iteration:
    k=k+1;
    
    % Berechnung des Fehlers in der 2-Norm fuer die Abbruchbedingung:
    fehler(k)=norm(f(x),2);
    
end

% Da k erhoeht wird, bevor der naechste Schritt tatsaechlich berechnet
% wurde, ziehen wir am Schluss wieder einen ab...:
k=k-1;