
% Lukáš Žúbor, 318061, 1.ročník, 2.stupeň, 3ZIJ1A 
clear;
clc;
syms x(n);
%system rovnic
e = [x(1)*[1;-1;-4],x(2)*[x(2);-x(2)^2;+4],x(3)*[-4;-4;1],[11;4;12]];
start_approx = [1;2;2];
precision = 3;
max_iterations = 100;
r = solve_nonlinear(e,precision,start_approx,max_iterations); %rovnica, presnost na desatinne miesta, pociatocne hodnoty, max pocet iteracii


function f = solve_nonlinear(B,precision,start_approx,max_iterations)
    syms x(n);
    syms j_helper;
    A = sym([]);
    disp("Solving system:");
    for i=1:length(B)-1
       A(i) = sum(B(i,:)); 
       disp(A(i));
    end
    size = length(A);
    jacobian = [];
    for i=1:size
        derived_line = [];
        for j=1:size
            A(i) = subs(A(i),x(j),j_helper);
            derived_line = [derived_line diff(A(i),x(j))];
            derived_line = subs(derived_line,j_helper,x(j));
        A(i) = subs(A(i),j_helper,x(j));
        end
        jacobian = [jacobian;derived_line];
    end
    for i=0:max_iterations
        disp("Iteration "+i);
        J = [];
        for j=1:size
            j_matrix_helper = [];
            for k=1:size
                j_matrix_helper = [j_matrix_helper subs(jacobian(j,k),x(k),start_approx(k))];
            end
        J = [J;j_matrix_helper];
        end
        F = [];
        for j=1:size
            f_matrix_helper = [];
            for k=1:size
                f_matrix_helper = [f_matrix_helper subs(B(j,k),x(k),start_approx(k))];
                if (k==size)
                   f_matrix_helper = [f_matrix_helper B(j,k+1)];
                end
            end
            F = [F;f_matrix_helper];
        end
        sum_f = [];
        for j=1:size
            sum_f = [sum_f; sum(F(j,:))];
        end
        sum_f = sum_f*-1;
        Newton = linsolve(J,sum_f);
        final = Newton+start_approx;
        if (round(final,precision)==round(start_approx,precision))
            disp("Done on "+num2str(i)+". iteration!");
            disp(round(vpa(start_approx),precision));
            f = vpa(start_approx);
            break
        end
        start_approx = final;
        f = vpa(start_approx);
        if (i == max_iterations)
            disp("I'm on the max iteration, but the result is not precise yet.");
            disp(vpa(start_approx));
        end
    end
end