% Solve the van der Pol equation and plot the solution in a
% parameter region where the equation is not very stiff.
%
% Thus, we can get away with using our fixed-step RK4.

global a;
T = 20;
N = 300;
a = 5;

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

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

plot (linspace(0,T,N+1),y(1,:), ";x;", ...
      linspace(0,T,N+1),y(2,:), ";y;"); 

