function yp = dwm(t,y,flag,delta,pp)
        mu = .01;                        % set parameter
%
        tint = ceil(t+.5./.5);
%
%       describe the system of differential equations
        yp1= mu.*y(1) - y(2).^2 +2.*y(3).^2 - delta.*y(3);
        yp2= y(2).*(y(1)-1)+pp(tint);
        yp3= mu.*y(3) +delta.*y(1) - 2.*y(1).*y(3);
%
%       evaluate function
        yp = [yp1;yp2;yp3];