What I am trying to do is solve a system of differential equations that involves a parameter δ0, then use that solution in another relation that modifies δ0. 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 = {Δ -> 10, γ -> 1, Ω -> 10, m -> 10^-25, ℏ -> 10^-34, k -> (2π)/10^-9};

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1;

{ρsol11, ρsol12, ρsol21, ρsol22} = 
    {ρ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, 0, 2}]; 

The goal is then to use the solution to ρsol22 in the following system:

fscatt = ℏ k^2 γ Re[ρsol22[t]]/m /. replace3;
δ0 = (2 π)/(500*10^-9)*10^3 - t fscatt;

Ultimately I would like to plot δ0 as a function of only t.

I imagined two strategies for doing this. The first was to actually solve the above system analytically. However, if you try running the above code you'll find that Mathematica is unable to execute the DSolveValue in a reasonable amount of time (or at all).

The second strategy was to assign δ0 an intial value and then change it incrementally according to the first equation in the above system using a fairly complicated Do loop. However, I ran into several problems with this approach, which you can read about here. Essentially I ran into a similar problem where Mathematica was taking forever.

Is there a better way to do this? I'm really at a loss at this point.


marked as duplicate by xzczd differential-equations Apr 15 at 15:07

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • $\begingroup$ Why not use ParametricNDSolve[]? $\endgroup$ – J. M. will be back soon Mar 31 at 22:22
  • $\begingroup$ @J.M.isaway would it be possible to do something like plot 'δ0' implicitly using 'ParametricNDSolve[]' or 'ParametricNDSolveValue[]' ? That might fix this issue. $\endgroup$ – Caleb Horwitz Apr 6 at 20:01
  • $\begingroup$ This is the same equation system as mentioned in your first question in this site, isn't it? Why don't you just use the analytic solution there? $\endgroup$ – xzczd Apr 10 at 3:31
  • $\begingroup$ And now the problem has been solved numerically, too. So this question has become a combination of the 2 linked questions, and I think it's time to mark it as duplicate. $\endgroup$ – xzczd Apr 15 at 15:11