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