% solve ODE system using RK4 with constant step size
function [x, derivativesTable] = rk4(functions, initialValues, interval, stepSize, maxSteps)
    x = initialValues;

    derivativesTable = buildDerivatiesTable(x, functions);
    
    % Calculate stepCount
    stepCount = ceil((interval(2) - interval(1)) / stepSize);
    if nargin == 5
        stepCount = min(stepCount, maxSteps - 1);
    end % IF we include max steps in our function input
    % (nargin is number of arguments in input)
    % then choose smaller number between maxSteps and stepCount and choose
    % it for stepCount

    [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, functions, derivativesTable);

    
    % append arguments to output
    x = [interval(1):stepSize:(stepCount * stepSize); x];
end

function derivativesTable = buildDerivatiesTable(x, functions)
    derivativesTable = zeros(size(x));
    for eqnum = 1:size(functions, 1)
        derivativesTable(eqnum, 1) = functions{eqnum}(x(:, 1));
    end
end

function [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, functions, derivativesTable)
    for step = 1 : stepCount
        [x, derivativesTable] = rk4stepLoop(x, step, functions, stepSize, derivativesTable);
    end
end

function [x, derivativesTable] = rk4stepLoop(x, step, functions, stepSize, derivativesTable)
    stepValue = x(:, step);
    
    for equationNumber = 1 : 2
        % generic single-step iteration
        phi = RK4Phi(functions{equationNumber}, stepValue, stepSize);
        x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * phi;
        
        % update derivatives table
        derivativesTable(equationNumber, step + 1) = functions{equationNumber}(x(:, step + 1));
    end
end