function result = aufgabe8()
    tolerance = 10^(-10);
    tolerance_mod = 10^(-12);
    x0 = [10;10;10];
    
    [x, steps, residuals_unmod, jacobi_residuals_unmod] = MyNewton(x0, tolerance, false);
    [x_mod, steps_mod, residuals_mod, jacobi_residuals_mod] = MyNewton(x0, tolerance_mod, true);
    
    x
    steps
    figure 1;
    semilogy(1:length(residuals_unmod), residuals_unmod);
    residuals_unmod
    figure 2;
    plot(1:length(jacobi_residuals_unmod), jacobi_residuals_unmod);
    
%    x_mod
%    steps_mod
%    figure 3;
%    semilogy(1:length(residuals_mod), residuals_mod);
%    residuals_mod
%    figure 4;
%    plot(1:length(jacobi_residuals_mod), jacobi_residuals_mod);
end

function [result, steps, residuals, jacobi_residuals] = MyNewton(x0, tolerance, modification)
    x_new = x0;
    x = [0;0;0];
    residual = 1;
    steps = 0;
    residual = 1;
    residuals = [];
    jacobi_residuals = [];
    
    while (residual > tolerance)
        x = x_new;
        [delta, jacobi_steps, jacobi_residuals_cur] = MyJacobi([0;0;0], DF(x(1), x(2), x(3)), -F(x(1), x(2), x(3)), tolerance, 5, modification);
%        delta = DF(x(1), x(2), x(3))\(-F(x(1), x(2), x(3)));
%        jacobi_steps = 0;
%        jacobi_residuals_cur = 0;
        x_new = x + delta;
        steps = steps + jacobi_steps;
        residual = norm(x_new-x);
        residuals = [residuals residual];
        jacobi_residuals = [jacobi_residuals jacobi_residuals_cur];
    end
    
    result = x_new;
end

function [delta, steps, residuals] = MyJacobi(x_ini, A, b, tolerance, iterations, modification)
    [n, m] = size(A);
    delta_new = x_ini;
    
    D = diag(diag(A));
    L = tril(A) - D;
    R = triu(A) - D;
    
    delta_cur = (zeros(n))(:,1);
    
    steps = 0;
    residual = 1;
    residuals = [];
    
    while (residual > tolerance)
        delta_cur = delta_new;
        delta_new = (D^(-1)) * (b - (L+R)*delta_cur);
        
        eig((D^(-1))*(D - A))
        pause
        
        steps = steps + 1;
        
        residual = norm(delta_new - delta_cur);
        residuals = [residuals residual];
        if (modification && steps >= 5)
            break
        end
    end
    
    delta = delta_new;
end

function zeroth = little_f(x,y,z)
    zeroth = 10*(x^2 + y^2 + z^2) + x*y^2*z + 2*x + y - (3/4)*z;
end

function first = F(x,y,z)
    first = [2 + 20*x + (y^2)*z; 1 + 20*y + 2*x*y*z; -(3/4) + x*y^2 + 20*z];
end

function second = DF(x,y,z)
    second = [20, 2*y*z, y^2; 2*y*z, 20 + 2*x*z, 2*x*y; y^2, 2*x*y, 20];
end