function [X] = RK4 (a, t0, T, N, f)
    h = T / N;
    p = size(a, 1);
    X = zeros(p, N);
    
    t = t0;
    X(:, 1) = a;
    k = h/6;
    l = h/2;
    
    for i=2:N
        K0 = f(t, X(:, i-1));
        K1 = f(t + l, X(:, i-1) + l * K0);
        K2 = f(t + l, X(:, i-1) + l * K1);
        K3 = f(t + h, X(:, i-1) + h * K2);
        X(:, i) = X(:, i-1) + k * ( K0 + 2 * K1 + 2 * K2 + K3 )
        t = t + h;
    end
    
endfunction