global a;
a = 5;
N = 10000;
T = 1000;

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

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

gset parametric
gsplot (y(:,2000:N))'