EDIT TO ADD: The below approach to depolymerization is wrong!
Hello MATLAB central,
I am trying to model a sequence of bond fission and radical recombination reactions in a polymer. In the below figure reaction (1) is the bond fission reaction and reaction (2) is the radical recombination reaction.
(Figure adapted from Levine et.al, 2009)
I am using kinetic equations of the form:
(1) where n is the length of the polymer A (of concentration [An]), [Rn] are the concentration of radicals R of length n, k1 and k2 are the reaction rates of the processes (2)
I am planning to use an ODE solver to get the concentration of all evolved species, but I am struggling with defining the sum of the products in reaction (1).
I first wanted to introduce additional matrix S that would describe products of all possible combinations of radicals that could produce the alkane of length j (see my code below), but it doesn't work. When using a stiff solver It returns a matrix of M with a warning:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.004156e-20.
While non-stiff solvers produce no warning message. The resulting matrix is just the initial conditions row repated for length of the timespan.
The code is: (please note that I used 23 as the length of the initial polymer, and I know I should probably have defined it as an additional argument of this function and I will change it after figuring out this issue).
function dMdt = levine(M,~,~, k1, k2, sumS, S)
dMdt = zeros(45,1);
for n = 1:45
for i =1:22
if n>=24 && n<=45 && i==1
S(i,n) = M(n).*M(69-n);
elseif n >=24 && n<=46-i && i>1 && i<=22
S(i,n) = M(n).*M(70-i-n);
sumS = sum(S);
if n>=1 && n<=22
dMdt(n) = -k1*M(n)+0.5*k2*sumS(n);
dMdt(n) = 0;
elseif n>23 && n<=45
dMdt(n) = -k1*M(n)+k2*sum(M(n-22:24));
And the command that I use to call this function is:
tspan = 1:1800;
m0 = [1 zeros(1,44)];
T = 673;
k1 = 1*10^16*exp(-89.7/(T*1.987*10^-3)) ;
k2 = 1.10*10^11*exp(-2.3/(T*1.987*10^-3)) ;
[t,M] = ode23s(@(t,M) levine(M,n,i,k1,k2,S,sumS), tspan, m0);