function [x,y] = RK(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N     
        k1 = f(x(k), y(k,:));
        k2 = f(x(k) + 3*h/4, y(k,:)+3*h*k1/4);
        
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1/3 + 2*k2/3)*h;        

    end;
end
