function yp=GWDf(t,y) 
global Pi alpha mu kappa beta gamma muW
yp(1,:)=Pi-beta.*y(1).*y(4)-mu.*y(1)+kappa.*y(3);
yp(2,:)=beta.*y(1).*y(4)-alpha.*y(2)-mu.*y(2);
yp(3,:)=alpha.*y(2)-kappa.*y(3)-mu.*y(3);
yp(4,:)=gamma.*y(3)-muW.*y(4);
yp(5,:)=gamma.*Pi./(mu.*(kappa+mu))-muW.*y(5);
