function [X,T] = RK4(x0,F,t,h)
% Implementation of the classical fourth order Runge-Kutta method for
% an autonomous system of ODEs in vector form
%
% INPUT: x0 (initial value vector)
%        F  (right hand side function)
%        t  (vector with initial and final time)
%        h  (step size)
%
% OUTPUT: X (matrix of approximate solutions at the discrete times)
%         T (vector of times t_i for which the solution is approximated)
%
n = round(diff(t)/h);
X = zeros(length(x0),n+1);
T = linspace(t(1),t(2),n+1);
X(:,1) = x0;
for i = 1:n
    K1 = h*F(X(:,i));
    K2 = h*F(X(:,i)+0.5*K1);
    K3 = h*F(X(:,i)+0.5*K2);
    K4 = h*F(X(:,i)+K3);
    X(:,i+1) = X(:,i) + (1/6)*(K1+2*K2+2*K3+K4);
end

end