function y=diff(t,u)

q = 0;
u_in =[0 0 0 0];
V = 100;

Et = u(1)+u(3);

% rate constants
k1   = 1/10;   
k_1  = 5.87/10; 
k2   = 0.69;  

% reaction rates
r1 = k1*u(1)*u(2); 
r2 = k_1*u(3);
r3 = k2*u(3); 

y=zeros(size(u));

y(1) = r2-r1+r3+q*(u_in(1)-u(1))/V;     %E
y(2) = r2-r1+q*(u_in(2)-u(2))/V;        %S
y(3) = r1-r2-r3+q*(u_in(3)-u(3))/V;     %ES
y(4) = r3+q*(u_in(4)-u(4))/V;           %P
y(5) = 0.6899*Et*u(2)/(12.7311+u(2)); %PM


%legend('E','S','ES','P')