
% Matúš Baláži, 318307, 1. ročník, 2. stupeň štúdia, 3ZIH1C

% parametre:
%   fun - sústava diferenciálnych rovníc
%   a - ľavá strana intervalu
%   b - pravá strana intervalu
%   h - veľkosť kroku
%   y_0 - vektor počiatočných podmienok pre y
function semestralna_praca_3_matus_balazi_318307(fun, a, b, h, y_0)

    if (nargin > 0 & nargin < 5)
        fprintf("Zadali ste nesprávny počet argumentov. " + ...
            "Zadajte všetky argumenty alebo ani jeden.\n");
        return;
    end

    % ak neboli zadané žiadne argumenty funkcie, 
    % tak sa použijú defaultné hodnoty
    if (nargin == 0)
        fun = @(t, y) [y(2); -y(1) - 2*exp(t) + 1; -y(1) - exp(t) + 1];
        a = 0;
        b = 1;
        h = 0.1;
        y_0 = [1; 0; 1];
    end

    t = a:h:b; % rozdelenie intervalu podľa kroku
    y = zeros(length(t), length(y_0)); % inicializácia vektora y so stĺpcom pre každú diferenciálnu rovnicu
    
    for i = 1:size(y_0)
        % nastavenie počiatočných podmienok pre y[1]
        y(1, i) = y_0(i);
    end

    k = zeros(4, numel(y_0)); % inicializácia vektorov k

    fprintf("<strong>\t\t%s\t\t%s\t\t\t\t\t\t%s</strong>\n", 'i', 't(i)', 'y(i)');
    
    for i = 1:(length(t))
        k(1, :) = fun(t(i), y(i,:));        
        k(2, :) = fun(t(i) + (h/2), y(i,:) + (h/2)*k(1,:));        
        k(3, :) = fun(t(i) + (h/2), y(i,:) + (h/2)*k(2,:));        
        k(4, :) = fun(t(i) + h, y(i,:) + h*k(3,:));
        
        y(i+1, :) = y(i, :) + (1 / 6) * (k(1, :) + (2 * k(2, :)) + (2 * k(3, :)) + k(4, :)) * h;
        
        fprintf('%9s %10s', num2str(i - 1), num2str(t(i)));
        for j = 1:length(y_0)
            fprintf('%15s', num2str(y(i, j)));
        end
        fprintf('\n')
    end

    fprintf('\nRiešenie:\n\t');
    for index = 1:length(y_0)
        fprintf('<strong>%7s</strong>    ', num2str(y(length(t), index)));
    end
    fprintf('\n');
end