# Help with complicated system of coupled differential equations with a parameter [duplicate]

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} =
DSolveValue[
{ρ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.

• Why not use ParametricNDSolve[]? Commented Mar 31, 2019 at 22:22
• @J.M.isaway would it be possible to do something like plot 'δ0' implicitly using 'ParametricNDSolve[]' or 'ParametricNDSolveValue[]' ? That might fix this issue.
– hclb
Commented Apr 6, 2019 at 20:01
• 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? Commented Apr 10, 2019 at 3:31
• 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. Commented Apr 15, 2019 at 15:11