Singular point problem in NDSolve

I am currently trying to solve a system of differential equations, but Mathematica returns an error. It seems like there is a problem with the third equation inside NDSolve. Is it because when L'[t] ==0, L''[t] cannot be calculated?

How can I solve this problem?

Qt = 10;
k1 = 83.5;
k2 = 0.0025;
A = 0.02^2 Pi;
l = 0.05;
T0 = 293;
R = 8.3145;
n = 0.02;
m = 0.05;
c = 0.493 m;
M = 0.05;


Above are constants

NDSolve[{Th'[t] == (Qt - (k1 A (Th[t] - Tc[t]))/l)/c,
Tc'[t] == ((k1 A (Th[t] - Tc[t]))/l - (k2 A (Tc[t] - T0))/L[t])/c,
M L''[t] L'[t] ==
Qt - D[5/2 n R (l Th[t] + (l + L[t]) Tc[t] + T0 L[t])/(
2 (l + L[t]))],
Th[0] == T0 + 20, Tc[0] == T0, L[0] == l, L'[0] == 8}, {Th, Tc,
L}, {t, 0, 1}]

• Welcome to Mathematica S.E. To start: 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, since the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) consider accepting the answer, if any, that solves your problem, by clicking checkmark sign, 4) give help too, by answering questions in your areas of expertise.
– bmf
Jan 27, 2023 at 0:05
• What error does it return? Jan 27, 2023 at 2:17
• It returns error "at t=~, step size is effectively zero; singularity or stiff system suspected." Jan 27, 2023 at 6:37

This might not be ideal solution but it gets rid of any warnings. By reducing the accuracy and increasing the number of steps and using "StiffnessSwitching" and this is important, need to Rationalize all parameters for this to work.

singularity or stiff system suspected are very hard to fix and it is case by case problem. Sometimes there is nothing you can do about singularity other than changing the mathematical model you used or changing the parameters used or initial conditions. From help:

Qt = 10 // Rationalize;
k1 = 83.5 // Rationalize;
k2 = 0.0025 // Rationalize;
A = 0.02^2 Pi // Rationalize;
l = 0.05 // Rationalize;
T0 = 293;
R = 8.3145 // Rationalize;
n = 0.02 // Rationalize;
m = 0.05 // Rationalize;
c = 0.493 m // Rationalize;
M = 0.05 // Rationalize;

ode1 = Th'[t] == (Qt - (k1 A (Th[t] - Tc[t]))/l)/c
ode2 = Tc'[t] == ((k1 A (Th[t] - Tc[t]))/l - (k2 A (Tc[t] - T0))/L[t])/c
ode3 = M L''[t] L'[t] == Qt - D[5/2 n R (l Th[t] + (l + L[t]) Tc[t] + T0 L[t])/(2 (l + L[t]))]

ic = {Th[0] == T0 + 20, Tc[0] == T0, L[0] == 1, L'[0] == 8}

sol = NDSolve[{ode1, ode2, ode3, ic}, {Th, Tc, L}, {t, 0, 1},
Method -> {"StiffnessSwitching"}, AccuracyGoal -> 5,
MaxSteps -> 10^7]


Plot[Evaluate[L[t] /. sol], {t, 0, 1}]


Plot[Evaluate[Th[t] /. sol], {t, 0, 1}]


It will be better not to use l as variable name as it looks like 1, and it is better to write each ode on its own, and write the initial conditions on its own. No need to throw everything into one command as that makes your code hard to read and modify.