function out=duffeq(t,xy,BwdOrFwd,vep,T)

% Francisco J. Beron-Vera, 2006/15/05

n=length(xy);
out=zeros(n,1);
x=xy(1:2:end-1);
y=xy(2:2:end);

out=[y;x-x.^3+vep*cos(2*pi*t/T)];
if strcmp(BwdOrFwd,'backward')
   out=-out;
end

out=reshape(out,[n/2 2])';
out=out(:);