Hi Hao,

You did a good job because solving ode is not easy job in Matlab. It took me a while to fix the bug in your code. I had to analyze your code with out plot and trying to display results after executing function in matlab, it lead me to error, which was declaration of function % Solve ODEs [T, Y] = ode45(@(t, Y) complexSystem(t, Y), tspan, Y0); The syntax supposed to be [t,y] = ode45(odefun,tspan,y0) For more information, please refer to https://www.mathworks.com/help/matlab/ref/ode45.html So, after updating this part of code, and using snippet code for plot, which still had error, plot(T, Y) should be replaced with plot(t,Y). Please try to execute code using modified code below and this should resolve your problem now. Also, for now use

% Initial conditions Y0 = [298; 0.75; 0.5; 0.5; 0.5; 0.5]; % Time span tspan = [0 10];

If you use 10000 for tspan, it will take longer for Matlab to execute and will never plot the function. Just another tip to keep in mind.

% Define the ODE function function dYdt = complexSystem(t, Y) % Constants and parameters Aa = 1e14; Ea = 2.24e-19; Ac = 4e13; Ec = 2.03e-19; As = 1e15; Es = 2.24e-19; Aec = 1e12; Eec = 1.4e-19; Kb = 1.38e-23; z0 = 0.033; 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);

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

% ODEs 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));

% Heat generation rates Qa = -ma * ha * dXadt; Qc = -mc * hca * dXcdt; Qs = -ma * hs * dXsdt; Qec = C * dSocdt; Q = Qa + Qc + Qs + Qec;

% Temperature change dTdt = Q;

% Output the derivatives dYdt = [dTdt; dXadt; dXcdt; dXsdt; dZdt; dSocdt]; end

% Initial conditions Y0 = [298; 0.75; 0.5; 0.5; 0.5; 0.5]; % Time span tspan = [0 10];

% Solve the ODEs [t, Y] = ode45(@complexSystem, tspan, Y0);

% Plot the results

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')

Please see attached snippet code.