
% Lukáš Žúbor, 318061, 1.ročník, 2.stupeň, 3ZIJ1A 
clear;
clc;
syms y(n);
%system rovnic, y(1) je x
e = y(1)^2-2*y(2)+y(3);
f = 2*y(1)^3+3*y(2)-y(3);
g = 2*y(1)^3+4*cos(y(2))-y(3);
initial_condition = [1;1;1];
interval = [0;0.1;0.5];
r = cauchy([e;f;g],initial_condition,interval); %rovnica, pociatocna podmienka, interval



function final = cauchy(equations,initial_condition,interval)
    syms y(t) ;
    step = interval(2);
    disp("Original_equations");
    disp(equations);
    
    for i = interval(1):interval(2):interval(3)-interval(2)
            disp("x(i):"+i);
            result = [];
            order = [];
            disp("Initial condition:");
            disp(initial_condition);
            %1. order
            old = [];
            for equations_count=1:length(equations)+1
                old = [old y(equations_count)];
            end
            new = [i];  
            for equations_count=1:length(equations)
                new = [new initial_condition(equations_count)];
            end
            for equations_count=1:length(equations)
                k(equations_count) = subs(equations(equations_count),old,new);
                order = [order k(equations_count)];
            end
            %2. order
            old = [];
            for equations_count=1:length(equations)+1
                old = [old y(equations_count)];
            end
            new = [i+step/2];  
            for equations_count=1:length(equations)
                new = [new initial_condition(equations_count)+step/2*order(length(order)-(length(equations)-equations_count))];
            end
            for equations_count=1:length(equations)
                k(equations_count) = subs(equations(equations_count),old,new);
                order = [order k(equations_count)];
            end
            %3. order
            old = [];
            for equations_count=1:length(equations)+1
                old = [old y(equations_count)];
            end
            new = [i+step/2];  
            for equations_count=1:length(equations)
                new = [new initial_condition(equations_count)+step/2*order(length(order)-(length(equations)-equations_count))];
            end
            for equations_count=1:length(equations)
                k(equations_count) = subs(equations(equations_count),old,new);
                order = [order k(equations_count)];
            end
            %4. order
            old = [];
            for equations_count=1:length(equations)+1
                old = [old y(equations_count)];
            end
            new = [i+step];  
            for equations_count=1:length(equations)
                new = [new initial_condition(equations_count)+step*order(length(order)-(length(equations)-equations_count))];
            end
            for equations_count=1:length(equations)
                k(equations_count) = subs(equations(equations_count),old,new);
                order = [order k(equations_count)];
            end
            for equations_count=1:length(equations)
                result = [result; initial_condition(equations_count) + (1/6)*step*(order(equations_count)+2*order(equations_count+length(equations))+2*order(equations_count+length(equations)+length(equations))+order(equations_count+length(equations)+length(equations)+length(equations)))];
            end
            initial_condition = vpa(result);
    end  
    final = vpa(result);
    disp("Final result in x = "+interval(3)+":");
    disp(round(final,3));

end