[T Y] = ode45('vdp1',0:1./60:20,[2; 0; 0; 0]);
figure;
plot(T,Y(:,1));