Coupled ODEs by ode45 and could not get the figures

Hi everyone!
I’ve been learning MATLAB for about six months, so I apologize if I ask any basic questions. In part of my project, I have a system of six coupled ODEs (all functions of t). You can refer to the attached picture for more details.
To obtain the T profile in terms of t, I decided to define the function similarly to how many others do:
function dYdt = complexSystem(t, Y)
Aa = 1e14; % Pre-exponential factor for xa
Ea = 2.24e-19; % Activation energy for xa
Ac = 4e13; % Pre-exponential factor for xc
Ec = 2.03e-19; % Activation energy for xc
As = 1e15; % Pre-exponential factor for xs
Es = 2.24e-19; % Activation energy for xs
Aec = 1e12; % Pre-exponential factor for xec
Eec = 1.4e-19; % Activation energy for xec
Kb = 1.38e-23; % Boltzmann constant
z0 = 0.033; % Initial value for z
ma = 0.13;
ha = 1.714e6;
mc = 0.29;
hca = 3.14e5;
hs = 2.57e5;
C = -3600 * 1.1 * 3.0 * (1-0.12-0.51);

T = Y(1);
Xa = Y(2);
Xc = Y(3);
Xs = Y(4);
Z = Y(5);
Soc = Y(6);

dXadt = -Aa * Xa * exp(-Ea/(Kb*T));
dXcdt = Ac * Xc * exp(1 – Xc) * exp(-Ec / (Kb * T));
dXsdt = -As * Xs * exp((-Z / z0) * (-Es / (Kb * T)));
dZdt = Aa * Xs * exp((-Z / z0) * (-Es / (Kb * T)));
dSocdt = -Aec * (1 – Xc) * Xa * exp(-Eec / (Kb * T));

Qa = -ma * ha * dXadt;
Qc = -mc * hca * dXcdt;
Qs = -ma * hs * dXsdt;
Qec = C * dSocdt;

Q = Qa + Qc + Qs + Qec;
dTdt = Q;

dYdt = [dTdt; dXadt; dXcdt; dXsdt; dZdt; dSocdt];
end
Then, I defined the initial conditions and solved the ODEs using ode45.
% Initial conditions
Y0 = [298; 0.75; 0.04; 0.015; 0.0331; 0]; % [T0; Xa0; Xc0; Xs0; Z0; Soc0]

% Time span
tspan = [0 100000];

% Solve ODEs
[T, Y] = ode45(@(t, Y) complexSystem(t, Y), tspan, Y0);
Lastly, I attempted to plot the results. However, the code took forever to run, and I didn’t get any figures from it.
Could my initial conditions be incorrect? Or are my ODEs not converging?
Thank you so much for your time and help!
figure;
plot(T, Y)
legend(‘T’, ‘X_a’, ‘X_c’, ‘X_s’, ‘Z’, ‘S_{oc}’)
title(‘Solution of Corrected Complex System’)
xlabel(‘Time t’)
ylabel(‘State Variables’)Hi everyone!
I’ve been learning MATLAB for about six months, so I apologize if I ask any basic questions. In part of my project, I have a system of six coupled ODEs (all functions of t). You can refer to the attached picture for more details.
To obtain the T profile in terms of t, I decided to define the function similarly to how many others do:
function dYdt = complexSystem(t, Y)
Aa = 1e14; % Pre-exponential factor for xa
Ea = 2.24e-19; % Activation energy for xa
Ac = 4e13; % Pre-exponential factor for xc
Ec = 2.03e-19; % Activation energy for xc
As = 1e15; % Pre-exponential factor for xs
Es = 2.24e-19; % Activation energy for xs
Aec = 1e12; % Pre-exponential factor for xec
Eec = 1.4e-19; % Activation energy for xec
Kb = 1.38e-23; % Boltzmann constant
z0 = 0.033; % Initial value for z
ma = 0.13;
ha = 1.714e6;
mc = 0.29;
hca = 3.14e5;
hs = 2.57e5;
C = -3600 * 1.1 * 3.0 * (1-0.12-0.51);

T = Y(1);
Xa = Y(2);
Xc = Y(3);
Xs = Y(4);
Z = Y(5);
Soc = Y(6);

Qa = -ma * ha * dXadt;
Qc = -mc * hca * dXcdt;
Qs = -ma * hs * dXsdt;
Qec = C * dSocdt;

Q = Qa + Qc + Qs + Qec;
dTdt = Q;

% Time span
tspan = [0 100000];

% Solve ODEs
[T, Y] = ode45(@(t, Y) complexSystem(t, Y), tspan, Y0);
Lastly, I attempted to plot the results. However, the code took forever to run, and I didn’t get any figures from it.
Could my initial conditions be incorrect? Or are my ODEs not converging?
Thank you so much for your time and help!
figure;
plot(T, Y)
legend(‘T’, ‘X_a’, ‘X_c’, ‘X_s’, ‘Z’, ‘S_{oc}’)
title(‘Solution of Corrected Complex System’)
xlabel(‘Time t’)
ylabel(‘State Variables’) Hi everyone!
I’ve been learning MATLAB for about six months, so I apologize if I ask any basic questions. In part of my project, I have a system of six coupled ODEs (all functions of t). You can refer to the attached picture for more details.
To obtain the T profile in terms of t, I decided to define the function similarly to how many others do:
function dYdt = complexSystem(t, Y)
Aa = 1e14; % Pre-exponential factor for xa
Ea = 2.24e-19; % Activation energy for xa
Ac = 4e13; % Pre-exponential factor for xc
Ec = 2.03e-19; % Activation energy for xc
As = 1e15; % Pre-exponential factor for xs
Es = 2.24e-19; % Activation energy for xs
Aec = 1e12; % Pre-exponential factor for xec
Eec = 1.4e-19; % Activation energy for xec
Kb = 1.38e-23; % Boltzmann constant
z0 = 0.033; % Initial value for z
ma = 0.13;
ha = 1.714e6;
mc = 0.29;
hca = 3.14e5;
hs = 2.57e5;
C = -3600 * 1.1 * 3.0 * (1-0.12-0.51);

T = Y(1);
Xa = Y(2);
Xc = Y(3);
Xs = Y(4);
Z = Y(5);
Soc = Y(6);

Qa = -ma * ha * dXadt;
Qc = -mc * hca * dXcdt;
Qs = -ma * hs * dXsdt;
Qec = C * dSocdt;

Q = Qa + Qc + Qs + Qec;
dTdt = Q;

% Time span
tspan = [0 100000];

Cart

Cart