function dy=seirf(t,y)
global beta sigma gamma f N

dy(1,1)=-beta*y(1)*y(3)/N;
dy(2,1)=beta*y(1)*y(3)/N-sigma*y(2);
dy(3,1)=sigma*y(2)-gamma*y(3);
dy(4,1)=(1-f)*gamma*y(3);