Free Radical Polymerization for Copolymers AA and Styrene using pseudokinetic rate constants

Hello, I have a MATLAB project for my polymer engineering class where I am supposed to use pseudokinetics to model the MW distribution and polymer concentration, etc of a styrene and acrylic acid copolymer. The first part of the project involved only acrylic acid and I will include the code below. Since this new case involves two polymers, I'm having a hard time and many sleepless nights making it into just one ode. I don't have much of a background in MATLAB (i have rarely used it prior to this semester, so ode45 is the only function I know for solving ODEs. If there's a better one I could use, I'd love to know, and I'd love any help with this project in general.

Part 1 code:

global kp kt kd f MWo

tf=15000;

tvect=[0 tf];

Mo=3.2; %mol/L changing monomer concentration +/- 25% changes polymer MW

Io=0.01; %mol/L

Mu=0;

MWo=72; %MWn=DPn*MWo, MWw=DPw*MWo

kp=950; %L/mol.s

Ep=6300; %cal/mol

T=348;%K

r=1.9872; %cal/K.mol

kt=110000000*exp(-2800/(r*T)); %L/mol.s

kd=140000000000000*exp(-30600/(r*T)); %s^-1

f=0.7;

%Ro=sqrt((2*f*kd*Io)/kt); %initial condition=Y1=Y2 ICs

%Y1= (kp*Mo/kt)+Yo; %initial condition

%Y2= ((2*kp*Mo*Y1)+(kp*Mo*Yo)+(kt*Yo^2))/(kt*Yo);

Ci0=[Mo Io Mu Mu Mu];

[t,Ci]=ode45(@monomer,tvect,Ci0);

M=Ci(:,1);I=Ci(:,2);muo=Ci(:,3);mu1=Ci(:,4);mu2=Ci(:,5);

%dmu1./dmu0=Y1./Yo;

%DPw=Y2./Y1;

%tau=kt.*Yo./(kp.*M);

Yo=sqrt((2*f*kd.*I)/kt);

tau=kt.*Yo./(kp.*M);

Y1= (kp.*M/kt)+Yo;

Y2=((2*kp.*M.*Y1)+(kp.*M.*Yo)+(kt.*Yo.^2))./(kt.*Yo);

MWn=(Y1./Yo)*MWo;

MWw=(Y2./Y1)*MWo;

DPnd=mu1./muo;

DPwd=mu2./mu1;

DPw=Y2./Y1;

DPn=Y1./Yo;

MWdn=DPnd*MWo;

MWdw=DPwd*MWo;

figure(1)

plot(t,Yo)

title('Yo vs t')

figure(6)

plot(t,Y1)

title('Y1 vs t')

figure(7)

plot(t,Y2)

title('Y2 vs t')

figure (10)

plot(t,muo)

title('Muo vs t')

figure (11)

plot(t,mu1)

title('Mu1 vs t')

figure (12)

plot(t,mu2)

title('Mu2 vs t')

figure(2)

plot(t,tau)

title('tau vs t')

figure(3)

plot(t,M)

title('[M] vs t')

figure(20)

plot(t,I)

title('[I] vs t')

figure(4)

plot(t,MWn,t,MWw)

legend('live MWn','live MWw')

title('live molecular weight distributions vs t')

figure(9)

plot(t,MWdn,t,MWdw)

legend('dead MWn','dead MWw')

title('dead molecular weight distributions vs time')

X=((Mo-M)./Mo)*100;

figure(15)

plot(t,X)

title('Monomer conversion vs time')

function dCdT=monomer(t,vars)

global kp kt kd f

M1=vars(1); I=vars(2);muo=vars(3); mu1=vars(4); mu2=vars(5);

dM=-kp*M*sqrt((2*f*kd.*I)/kt);

dI=-kd*I;

dmu0=2*f*kd*I;

dmu1=((kp*M/kt)+sqrt(2*f*kd.*I/kt))*kt.*sqrt((2*f*kd.*I)/kt);

dmu2=((2*kp.*M.*((kp.*M/kt)+sqrt(2*f*kd.*I/kt))+(kp.*M.*sqrt(2*f*kd.*I)/kt))+(kt.*(2*f*kd.*I/kt)))./(kt.*sqrt(2*f*kd.*I/kt))*kt.*sqrt(2*f*kd.*I/kt);

%muo=sqrt((2*f*kd.*I)/kt);

%mu1= (kp.*M/kt)+muo;

%mu2=((2*kp.*M.*mu1)+(kp.*M.*muo)+(kt.*muo.^2))./(kt.*muo);

%tau;

dCdT=[dM;dI;dmu0;dmu1;dmu2];

end

Part 2 (idk what I'm doing) code:

global kp11 kp12 kp21 kp22 ktc22 ktd11 ktd12 ktd22 f kd

tf=15000;

tvect=[0 tf];

M1o=1.6; %mol/L changing monomer concentration +/- 25% changes polymer MW

M2o=1.6; %styrene

Mo=M1o+M2o;

Io=0.01; %mol/L

Mu=0;

MW1=72; %MWn=DPn*MWo, MWw=DPw*MWo

MW2=104;

T = 348;

R = 1.987;

r1 = 0.27;

r2 = 0.72;

kp11 = 950;

kp12 = kp11/r1;

kp22 = (4.3e7)*exp(-7769/(R*T));

kp21 = kp22/r2;

ktc22 = (1.06e9)*exp(-1496/(R*T));

ktd11 = (1.1e8)*exp(-2800/(R*T));

ktd22 = (1.52e8)*exp(-1496/(R*T));

ktd12 = sqrt(ktd11*ktd22);

f = 0.7;

kd = (1.4e14)*exp(-30600/(R*T));

%Ro=sqrt((2*f*kd*Io)/kt); %initial condition=Y1=Y2 ICs

%Y1= (kp*Mo/kt)+Yo; %initial condition

%Y2= ((2*kp*Mo*Y1)+(kp*Mo*Yo)+(kt*Yo^2))/(kt*Yo);

Ci0=[Mo Io Mu Mu Mu];

[t,Ci]=ode45(@monomer,tvect,Ci0);

M=Ci(:,1);I=Ci(:,2);muo=Ci(:,3);mu1=Ci(:,4);mu2=Ci(:,5);

%dmu1./dmu0=Y1./Yo;

%DPw=Y2./Y1;

%tau=kt.*Yo./(kp.*M);

phi1 = (kp21.*M1)./((kp21.*M1)+(kp12.*f2.*M));

phi2 = 1-phi1;

f1 = M1./(M);

f2 = 1-f1;

M2=M-M1;

ktc = (ktc22.*phi2.^2);

ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);

Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));

kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);

Y1 = (kp.*M./(ktc+ktd))+Yo;

Y2 = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M/(ktc+ktd))+Yo;

tau=ktd.*Yo./(kp.*M);

beta=ktc*Yo./(kp.*M);

%M=M1+M2;

%MWn=(Y1./Yo)*MWo;

%MWw=(Y2./Y1)*MWo;

%DPnd=mu1./muo;

%DPwd=mu2./mu1;

%DPw=Y2./Y1;

%DPn=Y1./Yo;

%MWdn=DPnd*MWo;

%MWdw=DPwd*MWo;

%figure(1)

%plot(t,Yo)

%title('Yo vs t')

%figure(6)

%plot(t,Y1)

%title('Y1 vs t')

%figure(7)

%plot(t,Y2)

%title('Y2 vs t')

%figure (10)

%plot(t,muo)

%title('Muo vs t')

%figure (11)

%plot(t,mu1)

%title('Mu1 vs t')

%figure (12)

%plot(t,mu2)

%title('Mu2 vs t')

%figure(2)

%plot(t,tau)

%title('tau vs t')

%figure(3)

%plot(t,M)

%title('[M] vs t')

%figure(20)

%plot(t,I)

%title('[I] vs t')

%figure(4)

%plot(t,MWn,t,MWw)

%legend('live MWn','live MWw')

%title('live molecular weight distributions vs t')

%figure(9)

%plot(t,MWdn,t,MWdw)

%legend('dead MWn','dead MWw')

%title('dead molecular weight distributions vs time')

%X=((Mo-M)./Mo)*100;

%figure(15)

%plot(t,X)

%title('Monomer conversion vs time')

function dCdT=monomer(t,vars,kp11, kp12, kp21, kp22, ktc22, ktd11, ktd12, ktd22, f, kd)

M=vars(1);I=vars(2);%Y1=vars(8);Y2=vars(9);

dM=-((kp11*phi1*Yo*f1*M)+(kp12*phi1*Yo*f2*M)+(kp21*phi2*Yo*f1*M)+(kp22*phi2*Yo*f2*M));

dI=-kd*I;

dmu0=2*f*kd*I;

dmu1= kp*(M1+M2)*Y1*(tau+beta);

dmu2= kp*(M1+M2)*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));

%muo=sqrt((2*f*kd.*I)/kt);

%mu1= (kp.*M/kt)+muo;

%mu2=((2*kp.*M.*mu1)+(kp.*M.*muo)+(kt.*muo.^2))./(kt.*muo);

%tau;

dCdT=[dM;dI;dmu0;dmu1;dmu2];

end

[Davide Masiello]

One problem in your code is the way you handle passing extra-parameters to the ODE function.

Below, I have modified your part 1 code getting rid of global variables.

This is a much better approach.

Now, I suggest you to write a similar code for your part 2 and share it. Then, please explain a bit more extensively what changes in part 2 as compoared to part 1 and what is the problem.

Then we might try to come up with a solution.

Please remember to insert code as such, using the appropriate button in the menu above.

clear,clc

% Constants
Mo  = 3.2;                                          % mol/L changing monomer concentration +/- 25% changes polymer MW
Io  = 0.01;                                         % mol/L
Mu  = 0;
MWo = 72;                                           % MWn=DPn*MWo, MWw=DPw*MWo
kp  = 950;                                          % L/mol.s
Ep  = 6300;                                         % cal/mol
T   = 348;                                          % K
r   = 1.9872;                                       % cal/K.mol
kt  = 110000000*exp(-2800/(r*T));                   % L/mol.s
kd  = 140000000000000*exp(-30600/(r*T));            % s^-1
f   = 0.7;

% Numerical solution
Ci0     = [Mo Io Mu Mu Mu];
tvect   = [0 15000];
[t,Ci]  = ode45(@(t,vars)monomer(t,vars,kp,kt,kd,f),tvect,Ci0);

% Post-processing
M   = Ci(:,1);
I   = Ci(:,2);
muo = Ci(:,3);
mu1 = Ci(:,4);
mu2 = Ci(:,5);

Yo      = sqrt((2*f*kd.*I)/kt);
tau     = kt.*Yo./(kp.*M);
Y1      = (kp.*M/kt)+Yo;
Y2      = ((2*kp.*M.*Y1)+(kp.*M.*Yo)+(kt.*Yo.^2))./(kt.*Yo); 
MWn     = (Y1./Yo)*MWo;
MWw     = (Y2./Y1)*MWo;
DPnd    = mu1./muo;
DPwd    = mu2./mu1;
DPw     = Y2./Y1;
DPn     = Y1./Yo;
MWdn    = DPnd*MWo;
MWdw    = DPwd*MWo;
X       = ((Mo-M)./Mo)*100;

% Plots
figure(1)
subplot(3,4,1)
plot(t,Yo)
title('Yo vs t')
subplot(3,4,2)
plot(t,Y1)
title('Y1 vs t')
subplot(3,4,3)
plot(t,Y2)
title('Y2 vs t')
subplot(3,4,4)
plot(t,muo)
title('Muo vs t')
subplot(3,4,5)
plot(t,mu1)
title('Mu1 vs t')
subplot(3,4,6)
plot(t,mu2)
title('Mu2 vs t')
subplot(3,4,7)
plot(t,tau)
title('tau vs t')
subplot(3,4,8)
plot(t,M)
title('[M] vs t')
subplot(3,4,9)
plot(t,I)
title('[I] vs t')
subplot(3,4,10)
plot(t,MWn,t,MWw)
legend('live MWn','live MWw')
title('Live MW distributions')
subplot(3,4,11)
plot(t,MWdn,t,MWdw)
legend('dead MWn','dead MWw')
title('Dead MW distributions')
subplot(3,4,12)
plot(t,X)
title('Monomer conversion')

% ODE system
function dCdT=monomer(t,vars,kp,kt,kd,f)

M   =vars(1); 
I   =vars(2);
muo =vars(3); 
mu1 =vars(4); 
mu2 =vars(5);

dM   =-kp*M*sqrt((2*f*kd.*I)/kt);
dI   =-kd*I;
dmu0 =2*f*kd*I;
dmu1 =((kp*M/kt)+sqrt(2*f*kd.*I/kt))*kt.*sqrt((2*f*kd.*I)/kt);
dmu2 =((2*kp.*M.*((kp.*M/kt)+sqrt(2*f*kd.*I/kt))+(kp.*M.*sqrt(2*f*kd.*I)/kt))+(kt.*(2*f*kd.*I/kt)))./(kt.*sqrt(2*f*kd.*I/kt))*kt.*sqrt(2*f*kd.*I/kt);

dCdT = [dM;dI;dmu0;dmu1;dmu2];

end

[Oreoluwa Ogunnaike]

Thank you so much, I've edited my code to pass parameters using you suggestion.

A few changes in part 2: there are now two monomers, but just one dM ode, where M=M1+M2. there are also many new equations known as pseudokinetic constants (although they're not really constant at all), these ate the kp, ktc and ktd equations. dM= -kp*M*Yo. These equations are used in the ODE but they also depend on the ODE solution. My problem is, i dont know how to pass these variables into the ode, as well incorporating the M1 and M2 into the ode. I've included the code below. Thanks for helping out.

clear clc

%Constants
M1o=1.6; %mol/L changing monomer concentration +/- 25% changes polymer MW
M2o=1.6; %styrene
Mo=M1o+M2o;
Io=0.01; %mol/L
Mu=0;
MW1=72; %MWn=DPn*MWo, MWw=DPw*MWo
MW2=104;
T = 348; 
R = 1.987;
r1 = 0.27;
r2 = 0.72;
kp11 = 950;
kp12 = kp11/r1;
kp22 = (4.3e7)*exp(-7769/(R*T));
kp21 = kp22/r2;
ktc22 = (1.06e9)*exp(-1496/(R*T));
ktd11 = (1.1e8)*exp(-2800/(R*T));
ktd22 = (1.52e8)*exp(-1496/(R*T));
ktd12 = sqrt(ktd11*ktd22);
f = 0.7;
kd = (1.4e14)*exp(-30600/(R*T));

%Numerical solution 
Ci0=[Mo Io Mu Mu Mu];
tvect = [0 15000];
[t,Ci]=ode45(@(t,vars)monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f),tvect,Ci0);

% Post-processing
 M=Ci(:,1);
I=Ci(:,2);
muo=Ci(:,3);
mu1=Ci(:,4);
mu2=Ci(:,5);
%dmu1./dmu0=Y1./Yo;
%DPw=Y2./Y1;
%tau=kt.*Yo./(kp.*M);
phi1 = (kp21.*f1.*M)./((kp21.*f1.*M)+(kp12.*f2.*M));
phi2 = 1-phi1;
f1 = M1./(M);
f2 = 1-f1;
M2=M-M1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1 = (kp.*M./(ktc+ktd))+Yo;
Y2 = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M/(ktc+ktd))+Yo;
tau=ktd.*Yo./(kp.*M);
beta=ktc*Yo./(kp.*M);
%M=M1+M2;
%MWn=(Y1./Yo)*MWo;
%MWw=(Y2./Y1)*MWo;
%DPnd=mu1./muo;
%DPwd=mu2./mu1;
%DPw=Y2./Y1;
%DPn=Y1./Yo;
%MWdn=DPnd*MWo;
%MWdw=DPwd*MWo;
%figure(1)
%plot(t,Yo)
%title('Yo vs t')
%figure(6)
%plot(t,Y1)
%title('Y1 vs t')
%figure(7)
%plot(t,Y2)
%title('Y2 vs t')
%figure (10)
%plot(t,muo)
%title('Muo vs t')
%figure (11)
%plot(t,mu1)
%title('Mu1 vs t')
%figure (12)
%plot(t,mu2)
%title('Mu2 vs t')
%figure(2)
%plot(t,tau)
%title('tau vs t')
%figure(3)
%plot(t,M)
%title('[M] vs t')
%figure(20)
%plot(t,I)
%title('[I] vs t')
%figure(4)
%plot(t,MWn,t,MWw)
%legend('live MWn','live MWw')
%title('live molecular weight distributions vs t')
%figure(9)
%plot(t,MWdn,t,MWdw)
%legend('dead MWn','dead MWw')
%title('dead molecular weight distributions vs time')
%X=((Mo-M)./Mo)*100;
%figure(15)
%plot(t,X)
%title('Monomer conversion vs time')

%ODE system
function dCdT=monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f)
M=vars(1);I=vars(2);mu0=vars(3);mu1=vars(4);mu2=vars(5);
dM=-kp*M*Yo;
dI=-kd*I;
dmu0=2*f*kd*I;
dmu1= kp*M*Y1*(tau+beta);
dmu2= kp*M*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));
dCdT=[dM;dI;dmu0;dmu1;dmu2];

end

[Davide Masiello]

So, all the equations which depend on the integration variables are

phi1 = (kp21.*f1.*M)./((kp21.*f1.*M)+(kp12.*f2.*M));
phi2 = 1-phi1;
f1 = M1./(M);
f2 = 1-f1;
M2=M-M1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1 = (kp.*M./(ktc+ktd))+Yo;
Y2 = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M/(ktc+ktd))+Yo;
tau=ktd.*Yo./(kp.*M);
beta=ktc*Yo./(kp.*M);

Which you have listed under "post-processing" in your previous comment.

Am I correct?

[Oreoluwa Ogunnaike]

correct

[Davide Masiello]

I think we are close to writing a working code, but the problem is that by adding M1 and M2 as new variables you also need to add two new constraints (equations) in order for the problem to be solvable. One of these constraints is

M = M1+M2;

So there's still one to go.

[Oreoluwa Ogunnaike]

Hi Davide, I cannot find any other possible constraint, but it made me wonder why i couldnt have two separate ODEs for M1 and M2. I have adjusted the code for that. I still have the issue of figuring out how to use the equations in the ODE.

clear clc

%Constants
M1o=1.6; %mol/L changing monomer concentration +/- 25% changes polymer MW
M2o=1.6; %styrene
Io=0.01; %mol/L
Mu=0;
MW1=72; %MWn=DPn*MWo, MWw=DPw*MWo
MW2=104;
T = 348; 
R = 1.987;
r1 = 0.27;
r2 = 0.72;
kp11 = 950;
kp12 = kp11/r1;
kp22 = (4.3e7)*exp(-7769/(R*T));
kp21 = kp22/r2;
ktc22 = (1.06e9)*exp(-1496/(R*T));
ktd11 = (1.1e8)*exp(-2800/(R*T));
ktd22 = (1.52e8)*exp(-1496/(R*T));
ktd12 = sqrt(ktd11*ktd22);
f = 0.7;
kd = (1.4e14)*exp(-30600/(R*T));

%Numerical solution 
Ci0=[M1o M2o Io Mu Mu Mu];
tvect = [0 15000];
[t,Ci]=ode45(@(t,vars)monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f),tvect,Ci0);

% Post-processing
M1=Ci(:,1);
M2=Ci(:,2);
I=Ci(:,3);
muo=Ci(:,4);
mu1=Ci(:,5);
mu2=Ci(:,6);
%dmu1./dmu0=Y1./Yo;
%DPw=Y2./Y1;
%tau=kt.*Yo./(kp.*M);
phi1 = (kp21.*M1)./((kp21.*M1)+(kp12.*M2));
phi2 = 1-phi1;
f1 = M1./(M1+M2);
f2 = 1-f1;
ktc = (ktc22.*phi2.^2);
ktd = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1 = (kp.*(M1+M2)./(ktc+ktd))+Yo;
Y2 = (2.*kp.*(M1+M2).*Y1./((ktc+ktd).*Yo))+(kp.*(M1+M2)/(ktc+ktd))+Yo;
M=M1+M2;
tau=ktd.*Yo./(kp.*M);
beta=ktc*Yo./(kp.*M);
%MWn=(Y1./Yo)*MWo;
%MWw=(Y2./Y1)*MWo;
%DPnd=mu1./muo;
%DPwd=mu2./mu1;
%DPw=Y2./Y1;
%DPn=Y1./Yo;
%MWdn=DPnd*MWo;
%MWdw=DPwd*MWo;
%figure(1)
%plot(t,Yo)
%title('Yo vs t')
%figure(6)
%plot(t,Y1)
%title('Y1 vs t')
%figure(7)
%plot(t,Y2)
%title('Y2 vs t')
%figure (10)
%plot(t,muo)
%title('Muo vs t')
%figure (11)
%plot(t,mu1)
%title('Mu1 vs t')
%figure (12)
%plot(t,mu2)
%title('Mu2 vs t')
%figure(2)
%plot(t,tau)
%title('tau vs t')
%figure(3)
%plot(t,M)
%title('[M] vs t')
%figure(20)
%plot(t,I)
%title('[I] vs t')
%figure(4)
%plot(t,MWn,t,MWw)
%legend('live MWn','live MWw')
%title('live molecular weight distributions vs t')
%figure(9)
%plot(t,MWdn,t,MWdw)
%legend('dead MWn','dead MWw')
%title('dead molecular weight distributions vs time')
%X=((Mo-M)./Mo)*100;
%figure(15)
%plot(t,X)
%title('Monomer conversion vs time')

%ODE system
function dCdT=monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f)
M1=vars(1);M2=vars(2);I=vars(3);mu0=vars(4);mu1=vars(5);mu2=vars(6);
dM1=-kp*M1*Yo;
dM2=-kp*M2*Yo;
dI=-kd*I;
dmu0=2*f*kd*I;
dmu1= kp*(M1+M2)*Y1*(tau+beta);
dmu2= kp*(M1+M2)*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));
dCdT=[dM1;dM2;dI;dmu0;dmu1;dmu2];

end

[Davide Masiello]

Well, then it's easy enough.

I have attached an example in my answer below, so that you can click the accept button if you think it solves your problem.

Answer 1

[Davide Masiello]

You just need to calculate all the quantitites inside the ODE function as well.

If you need those quantities to be plotted, then you may recompute them after the integration as been performed (see the post-processing section in the example below).

clear, clc

% Constants
M1o     = 1.6; 
M2o     = 1.6; 
Mo      = M1o+M2o;
Io      = 0.01;
Mu      = 0;
MW1     = 72;
MW2     = 104;
T       = 348; 
R       = 1.987;
r1      = 0.27;
r2      = 0.72;
kp11    = 950;
kp12    = kp11/r1;
kp22    = (4.3e7)*exp(-7769/(R*T));
kp21    = kp22/r2;
ktc22   = (1.06e9)*exp(-1496/(R*T));
ktd11   = (1.1e8)*exp(-2800/(R*T));
ktd22   = (1.52e8)*exp(-1496/(R*T));
ktd12   = sqrt(ktd11*ktd22);
f       = 0.7;
kd      = (1.4e14)*exp(-30600/(R*T));

% Numerical solution 
Ci0     = [Mo Mo Io Mu Mu Mu];
tvect   = [0 15000];
[t,Ci]  = ode45(@(t,vars)monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f),tvect,Ci0);

% Post-processing
M1  = Ci(:,1);
M2  = Ci(:,2);
I   = Ci(:,3);
muo = Ci(:,4);
mu1 = Ci(:,5);
mu2 = Ci(:,6);

M       = M1+M2;
f1      = M1./(M);
f2      = 1-f1;
M2      = M-M1;

phi1    = (kp21.*f1.*M)./((kp21.*f1.*M)+(kp12.*f2.*M));
phi2    = 1-phi1;

ktc     = (ktc22.*phi2.^2);
ktd     = (ktd11.*phi1.^2)+(2*ktd12.*phi1.*phi2)+(ktd22.*phi2.^2);
Yo      = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp      = (kp11.*phi1.*f1)+(kp12.*phi1.*f2)+(kp21.*phi2.*f1)+(kp22.*phi2.*f2);
Y1      = (kp.*M./(ktc+ktd))+Yo;
Y2      = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M./(ktc+ktd))+Yo;
tau     = ktd.*Yo./(kp.*M);
beta    = ktc.*Yo./(kp.*M);

%ODE system
function dCdT=monomer(t,vars,kd,kp11,kp12,kp21,kp22,ktc22,ktd11,ktd12,ktd22,f)
M1  = vars(1);
M2  = vars(2);
I   = vars(3);
mu0 = vars(4);
mu1 = vars(5);
mu2 = vars(6);

M       = M1+M2;
f1      = M1./(M);
f2      = 1-f1;
M2      = M-M1;

phi1    = (kp21.*f1.*M)./((kp21.*f1.*M)+(kp12.*f2.*M));
phi2    = 1-phi1;

ktc     = (ktc22.*phi2.^2);
ktd     = (ktd11*phi1.^2)+(2*ktd12*phi1*phi2)+(ktd22*phi2.^2);
Yo      = sqrt((2.*f.*kd.*I)./(ktc+ktd));
kp      = (kp11*phi1*f1)+(kp12*phi1*f2)+(kp21*phi2*f1)+(kp22*phi2*f2);
Y1      = (kp.*M./(ktc+ktd))+Yo;
Y2      = (2.*kp.*M.*Y1./((ktc+ktd).*Yo))+(kp.*M/(ktc+ktd))+Yo;
tau     = ktd.*Yo./(kp.*M);
beta    = ktc*Yo./(kp.*M);


dM1   = -kp*M1*Yo;
dM2   = -kp*M2*Yo;
dI    = -kd*I;
dmu0  = 2*f*kd*I;
dmu1  = kp*M*Y1*(tau+beta);
dmu2  = kp*M*((tau*Y2)+(beta*(Y2+((Y1^2)/Yo))));
dCdT  = [dM1;dM2;dI;dmu0;dmu1;dmu2];
end

[Oreoluwa Ogunnaike]

This is perfect! Plotting graphs is easy enough for me. Thank you so much for all your help, I was so ready to give up.

[Davide Masiello]

Very happy to have helped