Several errors with InterpolatingFunction [closed]

Question

I'm getting some error messages that don't make sense. Apologies, as the code is somewhat complex. Basically what I'm trying to do is solve a system of differential equations, then fill a table with values for time and ρsol21[time], with a few other things going on. This is done via a Do loop.

h2 = {{0, -(Ω/2)}, {-(Ω/2), -k v + Δ}}

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

ρtderiv := -I (h2.ρ2 - ρ2.h2) + {{1/2 γ ρ22[t], -γ ρ12[t]}, {-γ ρ21[t], -(1/2) γ ρ22[t]}}

replace3 = {Δ -> 10, γ -> 1, Ω -> 10, m -> 10^-25, ℏ -> 1.0545718*10^-34, k -> (2 N[Pi, 10])/10^-9};

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

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = (2 N[Pi, 10])/10^-10;

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}]; 
 t0 += 0.01;
 δ0 -= 0.01*(2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3);
 txvarray[[i, 1]] = t0;
 txvarray[[i, 2]] = 2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3;
 txvarray[[i, 3]] = δ0;
 ρ120 = ρsol12[t0];
 ρ210 = ρsol21[t0];
 ρ220 = ρsol22[t0];
 ρ110 = ρsol11[t0],
 {i, 1, 200}]

I'm getting a few error messages that I don't quite understand. The first is

NDSolveValue::mxst: Maximum number of 10000 steps reached at the point t == 0.26002828028276154`.

And I don't have any idea what that means. The second is

InterpolatingFunction::dmval: Input value {0.27} lies outside the range of data in the interpolating function. Extrapolation will be used.

This doesn't make sense because when t0=0.27, the domain for the NDSolveValue of that iteration is {0.27,0.28}, so 0.27 should be in the range data. The third message is

An unknown box name (ToBoxes) was sent as the BoxForm for the expression. Check the format rules for the expression.

and I have no idea what this means either. Thanks for your time.

Answer 1

In principle, this is the working code. There are no mistakes, it is necessary only to define one parameter and avoid numerical divergence. All this is achieved if we put v=10^-9. Then we get a beautiful picture on the output

h2 = {{0, -(Ω/2)}, {-(Ω/2), -k*v + Δ}};

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

ρtderiv = 
  -I*(h2.ρ2 - ρ2.h2) + {{1/2*γ*ρ22[t], -γ*ρ12[t]}, {-γ*ρ21[t], -(1/2)*γ*ρ22[t]}};

replace3 = 
  {Δ -> 10, γ -> 1, Ω -> 10, m -> 10^-25, ℏ -> 1.0545718*10^-34, 
   k -> (2 N[Pi, 10])/10^-9, v -> 10^-9};

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

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = 2*Pi/10^-10;

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];
  t0 += 0.01;
  δ0 -= 0.01*(2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3);
  txvarray[[i, 1]] = t0;
  txvarray[[i, 2]] = 2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3;
  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, 2]]]}, {i, 1, 200}];
ListPlot[txabs, PlotStyle -> Orange]