%
% This function computes the vector field
% for the differential equations
%
%   x' = (x + y - x^3/3)/epsilon
%   y' = -x + a
%
function dvdt = vanderpol(t,v,p)
    %
    % Get the parameters from p
    % and the variables from v.
    %
    a       = p(1);
    epsilon = p(2);
    x = v(1);
    y = v(2);
   
    %
    % Compute the derivatives.
    %
    dxdt = (x + y - x^3/3)/epsilon;
    dydt = -x + a;
    
    %
    % Put the results into the vector
    % that is returned by this function.
    %
    dvdt = [dxdt; dydt];