function [ x, y, t ] = RK_2coupledODE( t0, tfin, N, x0, y0, f, g )

t=linspace(t0,tfin,N);
dt=t(2)-t(1);

x(1)=x0;
y(1)=y0;

% dx/dt=f(x,y);
% dy/dt=g(x,y);

for i=1:N-1
    xm(i)=x(i) + (dt/2)*( f(x(i), y(i)) );
    ym(i)=y(i) + (dt/2)*( g(x(i), y(i)) );
    
    x(i+1)=x(i) + dt*( f(xm(i), ym(i)) );
    y(i+1)=y(i) + dt*( g(xm(i), ym(i)) );
end

end

