global a;
T    = 7;
Nmax = 200;
a    = 0.21;
tol  = 10^-6;

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

[y,t] = ode_rk54 ('f', [0,T], [2;0], tol, Nmax);
yt = y';
save vanderpol1.dat yt

T     = 15;
[y,t] = ode_rk54 ('f', [0,T], [3;3.5], tol, Nmax);
yt = y';
save vanderpol2.dat yt

T     = 33;
[y,t] = ode_rk54 ('f', [0,T], [0.1;0.1], tol, Nmax);
yt = y';
save vanderpol3.dat yt

# plot (t, y(1,:), '-@', ...
#      t, y(2,:), '-@'); 

