Code Won't Finish running

Cant get my ODEs to run correctly with my X(2) value. Code will run but it doesn't complete and get a giant debugging page if I pause. Please help

[t,X] = ode45(@F,[0 340],[0; 400]);

figure
plot(t,X(:,1),'o')
legend('X')
xlabel('Reaction time (min)')
ylabel('X')

% Find t@X = 0.9
indx = find(abs(X(:,1)-0.9) < 0.01)
format
fprintf('Time (min) to achieve X = 0.9 is %f', t(indx))
fprintf('Conversion is %f', X(indx,1))

function dCdt = F(t,X)
R = 8.314; % J/mol/K
CA0 = 0.5; % M
CB0 = 0.75; % M
Cp = 3.8; % J/g/K
p = 0.75; %g/cm^3
dHrxn = -145000; % J/mol
V = 50; % L
NA0 = V*CA0; % mol
T0 = 400; %K
UA = 100; % W/K
ka = 1.4*10^7*exp(-7700/X(2)); % L/mol/min
dCdt = zeros(2,1);
dCdt(1,1) = ka*CA0*(1-X(1))*(1.5-X(1));
dCdt(2,1) = (-dHrxn/(p*Cp))*CA0^2*ka*(((1-X(1))*(1.5-X(1)))/1+(X(1)));
end

Answer 1

[Alan Stevens]

Use ode15s instead of ode45. Also, you might as well set the timespan to be [0 0.1] as nothing much changes after that!

Answer 2

[Mathieu NOE]

hello

your code works. I first reduced your time span because it was taking forever...

I also change the ode solver because - according to the plot - you shoud better use a solver for stiff ode.

speed up the whole thing quite significantly

the time transition you are looking for is around t = 0.02 s so no need to create a huge time span up to 340 s - unless there is something special happening at the end of the simulation ?

code more or less unchanged - minor stuff :

% [t,X] = ode45(@F,[0 340],[0; 400]);
[t,X] = ode15s(@F,[0 0.1],[0; 400]);

figure(1),
plot(t,X(:,1),'o')
legend('X')
xlabel('Reaction time (min)')
ylabel('X')
% Find t@X = 0.9
indx = find(abs(X(:,1)-0.9) < 0.01)
format
% fprintf('Time (min) to achieve X = 0.9 is %f', t(indx))
% fprintf('Conversion is %f', X(indx,1))

disp(['Time (min) to achieve X = 0.9 is :  ' num2str(t(indx)) '  s']);
disp(['Conversion is : ' num2str(X(indx,1)) ]);

function dCdt = F(t,X)
R = 8.314; % J/mol/K
CA0 = 0.5; % M
CB0 = 0.75; % M
Cp = 3.8; % J/g/K
p = 0.75; %g/cm^3
dHrxn = -145000; % J/mol
V = 50; % L
NA0 = V*CA0; % mol
T0 = 400; %K
UA = 100; % W/K
ka = 1.4*10^7*exp(-7700/X(2)); % L/mol/min

dCdt = zeros(2,1);
dCdt(1,1) = ka*CA0*(1-X(1))*(1.5-X(1));
dCdt(2,1) = (-dHrxn/(p*Cp))*CA0^2*ka*(((1-X(1))*(1.5-X(1)))/1+(X(1)));
% dCdt = [ka*CA0*(1-X(1))*(1.5-X(1));(-dHrxn/(p*Cp))*CA0^2*ka*(((1-X(1))*(1.5-X(1)))/1+(X(1)));];
%   % do the same as original code (the 3 line above are equivalent to this single line)

end