function [f, J] = numdiff(fun,x,delta)
% 
% 
%  [f, J] = numdiff(fun,x,delta) computes jacobian with difference quotient
%  
%  INPUT
%       fun     function handle
%       x       point x to evaluate f and approximate jacobian
%       delta   stepsize for difference approximation
%  OUTPUT
%       f       value f(x) 
%       J       jacobian at x
% 

% Übernehme Eingangsfunktion in f
f=fun(x);
% Gibt die Anzahl der Elemente im Array zurück und somit die Dimension von f.
m=numel(f);
% Gibt die Anzahl der Elemente im Array zurück und somit die Dimension von x.
n=numel(x); 

% Erstelle Jacobi-Matrix und fülle mit Nullen
J=zeros(m,n);
% Einheitsmatrix für e_j -> wird für Funktionsdefinition benötigt
I=eye(n);
% Schleife um die Jacobi-Matrix zu füllen.
for i=1:m
    for j=1:n
        % delta j berechnen
        dj     = delta(j);
        % Funktion definieren
        df     = (fun(x+dj*I(:,j))-fun(x));
        % Approximation berechnen und in Matrix speichern
        J(i,j) = df(i)/dj;
    end
end % function numdiff