function yp = vdp(t,y)
mu = 1;
yp = [y(2);
      mu*(1-y(1)^2)*y(2)-y(1)];