Cant get my loop to work

1 view (last 30 days)
Will Jeter
Will Jeter on 23 Nov 2020
Commented: Jan on 24 Nov 2020
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
  2 Comments
Mathieu NOE
Mathieu NOE on 24 Nov 2020
hello
no problem found
I run R2020b
Jan
Jan on 24 Nov 2020
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?

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 24 Nov 2020
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

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!