function [xdot,ydot]=fvdpol(t,x,y);

global a d e omega

% xdot=y;
% ydot=-x;
% r2=x^2+y^2;
% xdot= x/2-d*y/2 -x*r2/2;
% ydot=d*x/2+y/2-y*r2/2-a/4;
f=a*(1-omega^2)*cos(omega*t);

xdot=y;
ydot=-x-e*y*(x^2-1)+f;