# Plotting a function with a step increase from a delta function (top is cut off)

Using NDSolve I have two interpolating functions T[t] and n[t]. T[t] kind of includes a delta function to suddenly increase (will become clear below), but when plotting T[t] it has cut off the peak because it's too steep I think. How can I plot the full plot? Code as follows (starting with a load of jargon) Just defining a load of constants to be used in NDSolve:

{c1, c2, c3} = {4, 0.87, 0.5};
{kappa, zeta} = {1.7*10^-6, 1.2*10^-19};
L = 1*10^9;
kB = 1.38*10^-16;
a = 2*kappa/(35*kB*c2^(5/2)*c3*L^2);
b = (1 + c1)*zeta/(3*kB) - c1*c2*zeta/(5*c3*kB);
d = c1*c2*zeta/(5*c3*kB);
Lt = 1.2*10^70;
ti = -600;
tf = 4000;
T0 = 1*10^6;

n0 = (T0^2/(Lt*2*
c2^2))*(((((3*1.14*10^4)/7)^2 + 4*a*(c2^4)*(Lt^2)/d)^(1/
2)) - (3/7)*1.14*10^4);
R0 = (1.14*10^4)*(T0^2)/(n0*Lt*c2^2);
Qb = 3*kB*a*T0^(7/2)/(1 + (3/7)*R0) + 3*kB*b*n0^2/T0^(1/2);


Define the delta function:

qd = 7.5;
Qd = qd*DiracDelta[t - 25];
Qt = Qd + Qb;


Then I run the solver:

complete = Quiet@NDSolve[{T'[
t] == -(n[t]^-1) T[
t]^(7/2) (a) (1/(1 + (3/
7)*(((T[t]^2)/(c2^2*n[t]*(Lt)))*(1.14*10^4)))) -
n[t] T[t]^(-1/2) (b) + Qt/(3*kB*n[t]),
n'[t] ==
T[t]^(5/2) (a) (1/(1 + (3/
7)*((T[t]^2)/(c2^2*n[t]*(Lt))*(1.14*10^4)))) - (n[
t]^2) (T[t]^(-3/2)) (d), T[ti] == T0, n[ti] == n0}, {T,
n}, {t, ti, tf}, WorkingPrecision -> 40];


Now a) I know analytically that the peak of T[t] should be 1.51*10^7 and b) running FindMaximum returns a peak of 1.51*10^7. But when I plot the thing it only reaches about 1.4*10^7.

FindMaximum[T[t] /. complete, {t, 25, 26}, WorkingPrecision -> 20]

TPT = Plot[T[t] /. complete, {t, -50, 300}, PlotRange -> All,
Frame -> True, FrameLabel -> {"Time, s", "Temperature, K"}]


thanks!

## 3 Answers

In practice, what I would do is to inspect the peak more closely by doing a second plot:

TPT1 = Plot[T[t] /. complete, {t, 24, 26}, PlotRange -> All,
Frame -> True, FrameLabel -> {"Time, s", "Temperature, K"},
PlotRange -> All]


Then, having verified that this detail plot shows the correct peak height, I would combine it with the previous, less detailed plot:

Show[TPT, TPT1]


Since Show uses the options of the first plot in the argument list, you could also omit the Frame etc. from the plot TPT1.

{c1, c2, c3} = {4, 0.87, 0.5};
{kappa, zeta} = {1.7*10^-6, 1.2*10^-19};
L = 1*10^9;
kB = 1.38*10^-16;
a = 2*kappa/(35*kB*c2^(5/2)*c3*L^2);
b = (1 + c1)*zeta/(3*kB) - c1*c2*zeta/(5*c3*kB);
d = c1*c2*zeta/(5*c3*kB);
Lt = 1.2*10^70;
ti = -600;
tf = 4000;
T0 = 1*10^6;

n0 = (T0^2/(Lt*2*
c2^2))*(((((3*1.14*10^4)/7)^2 + 4*a*(c2^4)*(Lt^2)/d)^(1/2)) - (3/
7)*1.14*10^4);
R0 = (1.14*10^4)*(T0^2)/(n0*Lt*c2^2);
Qb = 3*kB*a*T0^(7/2)/(1 + (3/7)*R0) + 3*kB*b*n0^2/T0^(1/2);

qd = 7.5;
Qd = qd*DiracDelta[t - 25];
Qt = Qd + Qb;

complete =
Quiet@NDSolve[{T'[
t] == -(n[t]^-1) T[
t]^(7/2) (a) (1/(1 + (3/
7)*(((T[t]^2)/(c2^2*n[t]*(Lt)))*(1.14*10^4)))) -
n[t] T[t]^(-1/2) (b) + Qt/(3*kB*n[t]),
n'[t] == T[
t]^(5/2) (a) (1/(1 + (3/
7)*((T[t]^2)/(c2^2*n[t]*(Lt))*(1.14*10^4)))) - (n[t]^2) (T[
t]^(-3/2)) (d), T[ti] == T0, n[ti] == n0}, {T, n}, {t, ti, tf},
WorkingPrecision -> 40];

max = FindMaximum[T[t] /. complete, {t, 25, 26}, WorkingPrecision -> 20];

plateau = (T[t] /. complete /. t -> 25.)[[1]];


Rather than increasing the PlotPoints until the peak is displayed, just draw a line from the plateau to the max.

TPT = Plot[T[t] /. complete, {t, -50, 300},
PlotRange -> {Automatic, {0.5, 15.5} 10^6},
PlotStyle -> Blue,
Frame -> True,
Axes -> False,
FrameLabel -> {"Time, s", "Temperature, K"},
Epilog -> {Blue, Thick, Line[{{t, plateau}, {t, max[[1]]}}] /. max[[2]]}]


Add MaxRecursion->10 to Plot[...] :

{c1, c2, c3} = {4, 0.87, 0.5};
{kappa, zeta} = {1.7*10^-6, 1.2*10^-19};
L = 1*10^9;
kB = 1.38*10^-16;
a = 2*kappa/(35*kB*c2^(5/2)*c3*L^2);
b = (1 + c1)*zeta/(3*kB) - c1*c2*zeta/(5*c3*kB);
d = c1*c2*zeta/(5*c3*kB);
Lt = 1.2*10^70;
ti = -600;
tf = 4000;
T0 = 1*10^6;

n0 = (T0^2/(Lt*2*
c2^2))*(((((3*1.14*10^4)/7)^2 + 4*a*(c2^4)*(Lt^2)/d)^(1/
2)) - (3/7)*1.14*10^4);
R0 = (1.14*10^4)*(T0^2)/(n0*Lt*c2^2);
Qb = 3*kB*a*T0^(7/2)/(1 + (3/7)*R0) + 3*kB*b*n0^2/T0^(1/2);

qd = 7.5;
Qd = qd*DiracDelta[t - 25];
Qt = Qd + Qb;

complete = Quiet@NDSolve[{T'[
t] == -(n[t]^-1) T[
t]^(7/2) (a) (1/(1 + (3/
7)*(((T[t]^2)/(c2^2*n[t]*(Lt)))*(1.14*10^4)))) -
n[t] T[t]^(-1/2) (b) + Qt/(3*kB*n[t]),
n'[t] ==
T[t]^(5/2) (a) (1/(1 + (3/
7)*((T[t]^2)/(c2^2*n[t]*(Lt))*(1.14*10^4)))) - (n[
t]^2) (T[t]^(-3/2)) (d), T[ti] == T0, n[ti] == n0}, {T,
n}, {t, ti, tf}, WorkingPrecision -> 40];

TPT = Plot[T[t] /. complete, {t, -50, 300}, PlotRange -> All,
Frame -> True, FrameLabel -> {"Time, s", "Temperature, K"},MaxRecursion->10]


Here is the code to obtain the maximum seen by Plot[...,MaxRecursion->10] :

Cases[Normal[TPT],Line[___],{1,-1}] [[1,1,All]] //MaximalBy[#,Last]&


{{25.0047, 1.50174*10^7}}