function [x, y] = Heun(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
        %Predictor
        k1 = f(x(k),   y(k,:)       );
        k2 = f(x(k)+h, y(k,:) + k1*h);
        %corrector
        %disp([k1, k2, y + k1]);
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1 + k2 )*h/2.0;        
    end;
end

