function sys=md(x,u,K,flag)
% K P-control zero bias
oct=x(2); % states
holdup=x(1);
n=length(u)/2;
for i=1:n
    w(i)=u(2*(i-1)+1); %flow
    oct_inn(i)=u(2*i); % octane number
end
w=w(:);oct_inn=oct_inn(:);
wut=1/K*holdup;
if isequal(flag,'derivatives')
    sys(1)=sum(w)-wut;
    sys(2)=1/holdup*(w'*oct_inn - wut*oct -oct*sys(1)); 
elseif isequal(flag,'measurement')
    sys=[wut;oct;holdup];
else
    disp('ERROR')
    return
end
