function dy = SIV_eqn(t,y)
global B mu beta gamma alpha
dy = zeros(3,1); % need to be a column vector
dy(1) = B-mu*y(1) -beta *y(1)*y(3);
dy(2) = beta*y(1)*y(3)-mu*y(2);
dy(3) = gamma*y(2)-alpha*y(3);