function [x, y] = RK3(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) + h/2, y(k,:) + k1*h/2);
        k3 = f(x(k) + h,   y(k,:) - k1*h + 2*k2*h);
 
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1 + 4*k2 + k3)*h/6;        
    end;
end

