stiff ODE system plots disappear when I increase axes limits

Hello everyone! I recently asked for help to solve my ODE equations and I thought I might have to open another thread for this one.

My plots disappear when I set limits for them. I do see a line (straight line, probably due to the fact that it is only a tiny portion of the graph) before doing so, and nothing after adding limits on lines 24 and 28.

I need to be able to predict the rate of change until 5 min (300 second, with unit set on tspan as 0.001 s).

Thank you!

clc;
clear;
close all;

Q = @(v) sym(v);
%Units
cp = Q(29) ; %J/g*K
r = Q(8.314) ; %J/g*K
cv = cp - r;
tspan = [0 0.001] ; %sec
T0 = 300 ; %K
P0 = 20000000 ; %Pa (20 MPa in Pa due to R being in terms of m3)
syms T(t) P(t)
ode1 = diff(T) == 2.0*(T/P)*(r/cv)*(T0-T)-4*10^(-3)*(T/P)*(r^(2)/cv)*(P-P0);
ode2 = diff(P) == 2.0*(P/T)*(r/cv)*(T0-T)-4*10^(-3)*r*(cp/cv)*(P-P0);
odes = [ode1; ode2];
[eqs,vars] = reduceDifferentialOrder(odes, [T(t), P(t)]);
[M,F] = massMatrixForm(eqs,vars) ;
f = M\F ;
ode = odeFunction(f, vars) ;
[t,TP] = ode15s(ode, tspan, [T0+1, P0]);
plot(t, TP(:,1))
title('T')
axis([0 300000 120 300])
figure
plot(t, TP(:,2))
title('P')
axis([0 300000 0 20000000])
simplify(f)
subs(f,vars,[T0;P0])

Answer 1

[Abhiroop Rastogi]

Hi Michela,

The limits that you have chosen are very much out of bounds what the data is achieving. You can either let the limits be chosen by default or you can set the limits as shown below.

clc;
clear;
close all;
Q = @(v) sym(v);
%Units
cp = Q(29) ; %J/g*K
r = Q(8.314) ; %J/g*K
cv = cp - r;
tspan = [0 0.001] ; %sec
T0 = 300 ; %K
P0 = 20000000 ; %Pa (20 MPa in Pa due to R being in terms of m3)
syms T(t) P(t)
ode1 = diff(T) == 2.0*(T/P)*(r/cv)*(T0-T)-4*10^(-3)*(T/P)*(r^(2)/cv)*(P-P0);
ode2 = diff(P) == 2.0*(P/T)*(r/cv)*(T0-T)-4*10^(-3)*r*(cp/cv)*(P-P0);
odes = [ode1; ode2];
[eqs,vars] = reduceDifferentialOrder(odes, [T(t), P(t)]);
[M,F] = massMatrixForm(eqs,vars) ;
f = M\F ;
ode = odeFunction(f, vars) ;
[t,TP] = ode15s(ode, tspan, [T0+1, P0]);
plot(t, TP(:,1))
title('T')
%axis([0 300000 120 300])
axis([0 1e-3 300.999999992 301]) % modified limits

figure
plot(t, TP(:,2))
title('P')
%axis([0 300000 0 20000000])
axis([0 1e-3 1.999994e7 2e7]) % modified limits

simplify(f)

ans = 

subs(f,vars,[T0;P0])

ans =