global a;
global b;
a = 1;
b = 3;
N = 1000;
T = 100;

function y = f(t,x)
  global a;
  global b;
  y = [a - (b+1)*x(1) + x(1)^2 * x(2);
       b*x(1) - x(1)^2 * x(2)];
end

y = ode_rk4 ('f', [0,T], [0;0], N);

plot (y(1,:),y(2,:))