function yprime=duffn4(t,y,alpha,beta,delta,K,omega)
xx=y(1);
vv=y(2);
zz=y(3);
xsi=y(4);
ff1=vv;
ff2=-delta*vv-alpha*xx-beta*(xx^3)+K*zz;
ff3=xsi;
ff4=-(omega^2)*zz;
yprime=[ff1;ff2;ff3;ff4];