Cant get my loop to work

apparently there's something wrong in line 12 and i keep getting a debugging page for ode45 when i try to run it

% Initial conditions
Ca0 = 2;
Cb0 = 4;
Cc0 = 0;
Cd0 = 0;
T0  = 300:100:800;
IC = [repmat([Ca0 Cb0 Cc0 Cd0],length(T0),1),T0'];
% V range
Vspan = [0 10];
% Call ode solver
for k=1:size(IC,1)
 [V{k}, CT{k}] = ode45(@fn, Vspan, IC(k,:));
 figure;
 plot(V{k},CT{k}(:,1),V{k},CT{k}(:,2),V{k},CT{k}(:,3),V{k},CT{k}(:,4));
 xlabel('V (dm^3)'),ylabel('Ca.Cb.Cc,Cd (mol/dm^3)');
 legend('Ca','Cb','Cc','Cd');
 figure
 plot(V{k},CT{k}(:,5));
 grid on;
 xlabel('V(dm^3)'),ylabel('T(K)');
end
 
% ODE function
function dCTdV = fn(~,CT)
        % Data
        Cpa = 20;
        dh1a = 20000;
        dh2a = -10000;
        vo = 10;
        % Extract parameters
         Ca = CT(1);
         Cb = CT(2);
         Cc = CT(3);
         Cd = CT(4);
         T  = CT(5);
         
        % k and r functions
        k1a = 0.001*exp(-5000/1.987*(1/T - 1/300));
        k2a = 0.001*exp(-7500/1.987*(1/T - 1/300));
        r1a = -k1a*Ca*Cb^2;
        r2a = -k2a*Ca*Cc;
         % odes
         dCTdV = [(r1a+r2a)/vo;
                   2*r1a/vo;
                   (-2*r1a+r2a)/vo;
                   -2*r2a/vo;
                   (r1a*dh1a + r2a*dh2a)/((Ca+Cb+3*Cc+4*Cd)*vo*Cpa)];
end

[Mathieu NOE]

hello

no problem found

I run R2020b

[Jan]

I recommend to mention, why you assume, that there is a problem in line 12. Do you get an error message? Where does the debugger stop?

Answer 1

[Alan Stevens]

Replace your k loop by

for k=1:size(IC,1)
 [V, CT] = ode45(@fn, Vspan, IC(k,:));
 figure;
 plot(V,CT(:,1),V,CT(:,2),V,CT(:,3),V,CT(:,4));
 xlabel('V (dm^3)'),ylabel('Ca.Cb.Cc,Cd (mol/dm^3)');
 legend('Ca','Cb','Cc','Cd');
 figure
 plot(V,CT(:,5));
 grid on;
 xlabel('V(dm^3)'),ylabel('T(K)');
end