% the program takes the intial data x_0 and y_0 at the initial time t_0 and 
% computes up to time t_final using n timesteps.  It creates
% three vectors, t, x and y, which you can then plot against each other
% plot(x,y), hold on, plot(x(1),y(1),'o')
%
% [t,x,y] = sys_runge_kutta_2(x_0,y_0,t_0,t_final,n)

function [t,x,y] = sys_runge_kutta_2(x_0,y_0,t_0,t_final,n)

% size of the time-step:
dt = (t_final-t_0)/n;

% put initial data into vectors
t(1) = t_0;
x(1) = x_0;
y(1) = y_0;

for i=1:n

% careful, remember that f1 and f2 are _vectors_ and you need their
% first components for the dx/dt equation and their second components
% for the dy/dt equation.
  f1 = dt*F_sys(x(i),y(i),t(i));
  f2 = F_sys( x(i)+f1(1) , y(i)+f1(2), t(i) + dt); 

% do the time-stepping:
  x(i+1) = x(i) +  f1(1)/2 + dt*f2(1)/2;
  y(i+1) = y(i) +  f1(2)/2 + dt*f2(2)/2;
  t(i+1) = t_0 + i*dt;
end