To the best of my understanding you should include all three species in the ode-function. Perhaps something like this:
function xdot = f(t,x)
% Rate constants
k1 = 2.02;
% Reactants
cA = x(1);
cB = x(2);
cC = x(3);
% Reaction rates
r1 = k1*cA*cB;
% Differential equations
xdot = zeros(3,1);
xdot(1) = -r1;
xdot(2) = -r1;
xdot(3) = r1;
Then you specify the initial conditions for all 3 species:
%% Initial conditions
cA0 = 0.0109;
cB0 = 0.0109;
cC0 = 0;
x0 = [cA0; cB0; cC0];
%% Time span
tspan = [0 160]; % It seems necessary to increase time-span - slow reaction?
%% Solving ODE
[t,x] = ode45(@f,tspan,x0);
%% Plotting
ph = plot(t,x);
set(ph(1),'color','r');
set(ph(2),'color','b')
set(ph(3),'color','g','linewidth',2)
title('Ammonium Chloride Synthesis (Irreversible)');
xlabel('time');
Then I get curves that seem familiar to me - but I'm not an expert in this type of chemistry so I don't know if reaction-response should be on these time-scales or if they are completely off.
HTH
