Solving system of differential equations in loop with recursive variable

Question

Backslide introduced in 11.3, fixed in 12.0.

---

I've been posting various questions about similar pieces of code all year, so if this looks familiar you may have seen another one of my posts, but this is a unique issue that has not been asked yet, so I am posting it here. Basically I am using a loop so solve a system of coupled differential equations, using those results to change the value of a variable accordingly, then solve the system again, this time with the new value of the variable that was changed, and I want to do this in small steps 200 times, then plot the output. Here is the code:

h2 = {{0, -(Ω/2)}, {-(Ω/2), δ0 + Δ}};

ρ2 = {{ρ11[t], ρ12[t]}, {ρ21[t], ρ22[t]}};

ρdecay = {{1/2*γ*ρ22[t], -γ*ρ12[t]}, {-γ*ρ21[t], -(1/2)*γ*ρ22[t]}};

ρtderiv = -I*(h2.ρ2 - ρ2.h2) + ρdecay;

replace3 = {Δ -> -1*10^9, γ -> 1.6*10^9, Ω -> 1, m -> 10^-25, ℏ -> 1*10^-34, k -> (2 π)/(500*10^-9), v -> 10^3};

txvarray = Table[{0, 0, 0}, 200];

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = (2 π*10^3)/(500*10^-9);

Do[{ρsol11, ρsol12, ρsol21, ρsol22} = NDSolveValue[{ρ11'[t] == ρtderiv[[1, 1]], ρ12'[t] == ρtderiv[[1, 2]], ρ21'[t] == ρtderiv[[2, 1]], ρ22'[t] == ρtderiv[[2, 2]], ρ11[t0] == ρ110, ρ12[t0] == ρ120, ρ21[t0] == ρ210, ρ22[t0] == ρ220} /. replace3, {ρ11, ρ12, ρ21, ρ22}, {t, t0, t0 + 0.01}, MaxSteps -> Infinity];
Δt = 0.01;
t0 += Δt;
fscatt = ℏ k^2 γ Re[ρsol22[t0]]/m /. replace3;
δ0 -= Δt fscatt;
txvarray[[i, 1]] = t0;
txvarray[[i, 2]] = fscatt;
txvarray[[i, 3]] = δ0;
ρ120 = ρsol12[t0];
ρ210 = ρsol21[t0];
ρ220 = ρsol22[t0];
ρ110 = ρsol11[t0], {i, 1, 200}]

txabs = Table[{txvarray[[i, 1]], Abs[txvarray[[i, 3]]]}, {i, 1, 200}];

ListPlot[txabs, PlotStyle -> Green]

Note the δ0 in the first line, in h2. This is the variable that changes after each iteration (this happens in the line δ0 -= Δt fscatt;). The problem is that whenever I try to run this code, it gets stuck and I have to quit the kernal. Any ideas why this is happening/how to fix it? Thanks!

Answer 1

Update

Just tested on Wolfram cloud, the sample no longer gets stuck in v12.0.

---

OK, I manage to reproduce the issue in v11.3, but not in v11.2 and v9. I guess this is related to the truth that "CatchMachineUnderflow" option is removed in this version, because after adding WorkingPrecision -> 16 to NDSolveValue together with Rationalize[..., 0], the problem is resolved:

h2 = {{0, -(Ω/2)}, {-(Ω/2), δ0 + Δ}};   
ρ2 = {{ρ11[t], ρ12[t]}, {ρ21[t], ρ22[t]}};
ρdecay = {{1/2γρ22[t], -γρ12[t]}, {-γρ21[t], -(1/2)γρ22[t]}};
ρtderiv = -I*(h2.ρ2 - ρ2.h2) + ρdecay;
replace3 = {Δ -> -110^9, γ -> 1.610^9, Ω -> 1, 
   m -> 10^-25, ℏ -> 110^-34, k -> (2 π)/(50010^-9), v -> 10^3};
txvarray = Table[{0, 0, 0}, 200];
t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = (2 π10^3)/(50010^-9);
Do[{ρsol11, ρsol12, ρsol21, ρsol22} = 
  NDSolveValue[{ρ11'[t] == ρtderiv[[1, 1]], ρ12'[
        t] == ρtderiv[[1, 2]], ρ21'[t] == ρtderiv[[2, 1]], ρ22'[
        t] == ρtderiv[[2, 2]], ρ11[t0] == ρ110, ρ12[
        t0] == ρ120, ρ21[t0] == ρ210, ρ22[t0] == ρ220} /. 
     replace3 // Rationalize[#, 0] &, {ρ11, ρ12, ρ21, ρ22}, {t, t0, 
    t0 + 0.01}, WorkingPrecision -> 16];
 Δt = 0.01;
 t0 += Δt;
 fscatt = ℏ k^2 γ Re[ρsol22[t0]]/m /. replace3;
 δ0 -= Δt fscatt;
 txvarray[[i, 1]] = t0;
 txvarray[[i, 2]] = fscatt;
 txvarray[[i, 3]] = δ0;
 ρ120 = ρsol12[t0];
 ρ210 = ρsol21[t0];
 ρ220 = ρsol22[t0];
 ρ110 = ρsol11[t0], {i, 1, 200}]
txabs = Table[{txvarray[[i, 1]], Abs[txvarray[[i, 3]]]}, {i, 1, 200}];
ListPlot[txabs, PlotStyle -> Green]

The output is just the same as in v9 and v11.2:

v11.3 result after adding WorkingPrecision -> 16:

v9 result (Definition of txvarray is modified to txvarray = Table[{0, 0, 0}, {200}]):