function dc = fzad2(~,c)

global k1 k2 k3 k4 k5 k6 k7 k8 k9

dc=zeros(9,1);

k1=2.95e-1; %[1/(M*s)]%
k2=2e10;    %[1/(M*s)]%
k3=1e5;     %[1/(M*s)]%
k4=3e5;     %[1/(M^3*s)]%
k5=3e12;    %[1/(M^2*s)]%
k6=2.2;     %[1/s]%
k7=1e6;     %[1/(M*s)]%
k8=5e10;    %[1/(M*s)]%
k9=3e3;     %[1/(s)]%

%C(1): [IO3-] C(2): [HSO3-] C(3): [HIO2] C(4): [I-] C(5): [H+] C(6): [HOI]%
%c(7): [I2] c(8): [SO3(2-)] c(9): [SO4(2-)] %

r(1)=k1*c(1)*c(2);             % IO3(-)+HSO3(-)→SO4(2-)+HIO2 %
r(2)=k2*c(3)*c(4)*c(5);        % HIO2+I(-)+H(+)→2HOI %
r(3)=k3*c(3)*c(6);             % HIO2+HOI→IO3(-)+I(-)+2H(+) %
r(4)=k4*c(1)*c(4)*(c(5))^2;    % IO3(-)+I(-)+2H(+)→HIO2+HOI %
r(5)=k5*c(6)*c(4)*c(5);        % HOI+I(-)+H(+)→I2+H2O %
r(6)=k6*c(7)*c(8);             % I2+H2O→HOI+I(-)+H(+) %
r(7)=k7*c(7)*c(2)*c(8);        % I2+HSO3(-)+H2O→2I(-)+SO4(2-)+3H(+) %
r(8)=k8*c(9)*c(5);             % SO3(2-)+H(+)→HSO3(-) %
r(9)=k9*c(2);                  % HSO3(-)→SO3(2-)+H(+) %

dc(1)=-r(1)+r(3)-r(4);                                  % [IO3-] %
dc(2)=-r(1)-r(7)+r(8)-r(9);                             % [HSO3-] %
dc(3)=-r(2)-r(1)+r(4);                                  % [HIO2] %
dc(4)=-r(2)+r(3)-r(4)-r(5)+r(6)+2*r(7);                 % [I-] %
dc(5)=-r(2)+2*r(3)-2*r(4)-r(5)+r(6)+3*r(7)-r(8)+r(9);   % [H+] %
dc(6)=2*r(2)-r(3)+r(4)-r(5)+r(6);                       % [HOI] %
dc(7)=r(5)-r(6)-r(7);                                   % [I2] %
dc(8)=-r(8)+r(9);                                       % [SO3(2-)] %
dc(9)=r(1);                                             % [SO4(2-)] %
end




clear all
clc

global k1 k2 k3 k4 k5 k6 k7 k8 k9
%C(1): [IO3-] C(2): [HSO3-] C(3): [HIO2] C(4): [I-] C(5): [H+] C(6): [HOI]%
%c(7): [I2] c(8): [SO3(2-)] c(9): [SO4(2-)] %

c0 = [0.0715; 0.0101; 0.0; 1e-6; 8.36e-9; 0.0; 0.0; 0.0764; 0.0];
tspan= [0,480]; % [s] %
[t, c]=ode45(@(t,c)fzad2 ,tspan,c0)