function y=ozonea(u)
% Decomposition of ozone. 
% Atkins and Jones, first edition, p. 626
% This function computes the quotient 
% between the rate u3' and the overall rate
% law -u3^2/u2. The quotient should be 
% 2*k21*k11/k12.

k11=10;
k12=10;
k21=.01;
% k22=.0001;
k22=0;

% k11=1;k12=1;k21=1;

r11=k11*u(:,3);
r12=k12*u(:,2).*u(:,1);
r21=k21*u(:,1).*u(:,3);
r22=k22*u(:,2).^2;

% y=-(-r11+r12-r21+r22) ./ ( (u(:,3).^2) ./ u(:,2) );

y(:,1)=( -(2*k21*k11)*( (u(:,3).^2)+eps) ./ ( k12*u(:,2)+k21*u(:,3) +eps) )-(-r11+r12-r21+r22+eps);

y(:,2)=( -(2*k21*k11)*( (u(:,3).^2)+eps) ./ ( k12*u(:,2)+k21*u(:,3) +eps) );

y(:,3)=(-r11+r12-r21+r22+eps);

y(:,4)=( -(2*k21*k11)*( (u(:,3).^2)+eps) ./ ( k12*u(:,2)+k21*u(:,3) +eps) )./(-r11+r12-r21+r22+eps);
