function y=diff(t,u)

% rate constants 
K1  = 1.2e-2;   %R2
k12 = 1e10;
K2  = 7.4e9;    %R3
k21 = 1.4e11;

% reaction rates 
r11 = k12*K1*u(1);
r12 = k12*u(2)*u(3);
r21 = k21*u(2)*u(4);
r22 = k21/K2*u(5);

y=zeros(size(u));

y(1) = r12-r11;           %HSO4  
y(2) = r11-r12+r22-r21;   %H+ 
y(3) = r11-r12;           %SO4
y(4) = r22-r21;           %OH-
y(5) = r21-r22;           %H20

%legend('HSO_4^-','H^+','SO_4^2^-','OH^-','H_2O')
%plot(t,-log10(u(:,2)))