ODE45has long runtime and graph will not plot

[Walter Roberson]

It is probably a stiff system; use ode15s instead of ode45

Answer 1

[Alan Stevens]

Here's my attempt to make sense of your equations

tauspan = 0:100:10000;

%  p0 = [0.5 0.5 0 0];
%  dp1dt = -rf + rr
%  dp2dt = 2*(-rf + rr)
%  dp3dt = rf - rr
%  dp4dt = 4*(rf - rr)
%  d(2*p1-p2)/dt = 0   so  2*p1 - p2 = c2   (where c2 is a constant)
%  d(p1 + p3)/dt = 0   so    p1 + p3 = c3   (where c3 is a constant)
%  d(4*p1+p4)/dt = 0   so  4*p1 + p4 = c4   (where c4 is a constant)
% Using initial conditions we have:
% 2*0.5 - 0.5 = c2     so  p2 = 2*p1 - 0.5
% 0.5 + 0 = c3         so  p3 = -p1 + 0.5
% 4*0.5 + 0 = c4       so  p4 = -4*p1 + 2

kf = 1.32E10*exp(-236.7)/8.314;
% kr = 1.32E10*exp(-285.2)/8.314
% rf = kf*p1*(2*p1-0.5)^2/T
% rr = kr*(-p1+0.5)*(-4*p1+2)^2/T

% Scale the time base:
% tau = sf*t    dp1dt = dp1dtau*dtau/dt = dp1dtau*sf
% sf*dp1dtau = -kf*p1*(2*p1-0.5)^2/T + kr*(-p1+0.5)*(-4*p1+2)^2/T
% Let kf/sf = 1  so 
sf = kf;
% and  kr/sf = exp(-285.2)/exp(-236.7) = exp(-48.5)
% dp1dtau = (-p1*(2*p1-0.5)^2 + exp(-48.5)*(-p1+0.5)*(-4*p1+2)^2)/T

p10 = 0.5;
T = 900:75:1200;

for i = 1:numel(T)
[tau, p1] = ode45(@(t,p) ODET(t,p,T(i)), tauspan, p10);
p2 = 2*p1-0.5;
p3 = -p1+0.5;
p4 = -4*p1+2;
t = tau/sf;
figure
hold on
Tlbl = ['T = ', int2str(T(i))];
subplot(2,2,1)
plot(t,p1),grid
title(Tlbl)
xlabel('t'), ylabel('p1')
subplot(2,2,2)
plot(t,p2),grid
title(Tlbl)
xlabel('t'), ylabel('p2')
subplot(2,2,3)
plot(t,p3),grid
title(Tlbl)
xlabel('t'), ylabel('p3')
subplot(2,2,4)
plot(t,p4),grid
title(Tlbl)
xlabel('t'), ylabel('p4')
hold off
end

function dpdtau = ODET(~,p,T)

         dpdtau = (-p*(2*p-0.5)^2 + exp(-48.5)*(-p+0.5)*(-4*p+2)^4)/T;
        
end