function dy = hum_mos_eqn(t,y)
global B Bp mu beta mup mubar
dy = zeros(4,1); % need to be a column vector
dy(1) = B-mu*y(1) -beta *y(1)*y(4);
dy(2) = beta*y(1)*y(4)-mubar*y(2);
dy(3) = Bp-mup*y(3)-beta*y(3)*y(2);
dy(4) = beta*y(3)*y(2)-mup*y(4);