MATLAB: Problem using ODE15s

chemistryequilibrium reactionkineticsodereactions in series

Hey everyone. I'm trying to model this chemical reaction:
This is a combination of equilibrium reactions (forward and backward reactions) and reactions in series.
I have written a code that does not work really well, the CO2 concentration seems ok but the other concentrations are just step functions and it is not supposed to look like that. My main guess is that something is wrong with the rates of reaction but I cannot see what is wrong. My second guess is that I have chosen wrong coding strategy.
clc
close all
clear all
c_0 = (46./1e6).*1000; %mol/m^3



H2CO3_0=0; %mol/m^3
HCO3_0=6437e-6*1000; %mol/m^3
H_0=0.00000001*1000; %mol/m^3
tspan= [0 1000];
u0=[c_0;H2CO3_0;HCO3_0;H_0]; %initial conditions for ODE's
%diff equations----------------------------
%u(1)=concentration CO2
%u(2) concentration H2CO3
%u(3) concentration HCO3-
%u(4) concentration H+
k_a=18; %s^-1 Backward constant
k_b=0.04; %s^-1 Forward constant

k_21=1e7; %s^-1 Forward constant
k_12=5e10/1000; %m^3/(mol*s) Backward constant
pak1=@(t,u) [(k_a.*u(2))-(k_b.*u(1)); %dCO2/dt
    ((k_12.*u(3).*u(4))-k_21.*u(2)); %dH2CO3/dt
    -((k_12.*u(3).*u(4))-k_21.*u(2)); %dHCO3/dt
    -((k_12.*u(3).*u(4))-k_21.*u(2))]; %dH/dt
opt=odeset('AbsTol',1e-10,'RelTol',1e-12);
[t u] = ode15s(pak1, [tspan], u0, opt);
subplot(2,2,1)
plot(t./3600,u(:,1).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('CO$_2$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
subplot(2,2,2)
plot(t./3600,u(:,2).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('H$_2$CO$_3$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
subplot(2,2,3)
plot(t./3600,u(:,3).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('HCO$_3^{-1}$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
subplot(2,2,4)
plot(t./3600,-log10(u(:,4)./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('pH','Interpreter','LaTeX');

Best Answer

You have an error in your continuity-equation, you have forgotten the loss of H2CO3 due to the first backwards reaction H2CO3 -> CO2 + H20: k_a, and you haven't included the source due to the CO2 + H2O: k_b either. I think it should look like this:
pak1=@(t,u) [(k_a.*u(2))-(k_b.*u(1)); %dCO2/dt
    ((k_12.*u(3).*u(4))+k_b.*u(1)-(k_21 + k_a).*u(2)); %dH2CO3/dt
    -((k_12.*u(3).*u(4))-k_21.*u(2)); %dHCO3/dt
    -((k_12.*u(3).*u(4))-k_21.*u(2))]; %dH/dt
I haven't checked for more. Your problem is painfully familiar to me, in order to get these right I typically have to print out the code, and explicitly mark the source and loss-terms of each reaction making sure that they are properly included in the reactions for each reactant and product (multicoloured pens help.)
HTH
Related Question