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
% 

f = fun(x); %Funktionswert berechnen
I = eye(size(x,1)); %einheitsmatrix für ej
J = zeros(size(f,1),size(x,1)); %Jacobimatrix in richtiger Größe initialisieren

for i=1:size(f)
    for j=1:size(x)
        deltaj = delta(j); %Aktuelle Schrittweite holen
        deltaf = fun(x+deltaj*I(:,j))-fun(x); %Funktion berechnen
        J(i,j) = deltaf(i)/deltaj; %Matrixkomponente berechnen und speichern
    end
end

end % function numdiff