1
$\begingroup$

I want to solve the following time dependent coupled differential equations numerically. Please guide me.

{
 {Derivative[1][Subscript[\[Rho], 1, 1]][t] == \[Gamma]/
    2 - \[Gamma] Subscript[\[Rho], 1, 1][t] - 
    I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] - 
       1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t]) + 
    1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t]},
 {Derivative[1][Subscript[\[Rho], 1, 2]][
    t] == -\[Gamma] Subscript[\[Rho], 1, 2][t] - 
    I (\[CapitalDelta]a Subscript[\[Rho], 1, 2][
         t] - \[CapitalDelta]b Subscript[\[Rho], 1, 2][t] + 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 3][t] - 
       1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 2][t])},
 {Derivative[1][Subscript[\[Rho], 1, 3]][t] == 
   1/2 (-\[Gamma] Subscript[\[Rho], 1, 3][
         t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 1, 3][
         t]) - I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] + 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 2][
         t] + \[CapitalDelta]a Subscript[\[Rho], 1, 3][t] - 
       1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t])},
 {Derivative[1][Subscript[\[Rho], 2, 1]][
    t] == -\[Gamma] Subscript[\[Rho], 2, 1][t] - 
    I (-\[CapitalDelta]a Subscript[\[Rho], 2, 1][
         t] + \[CapitalDelta]b Subscript[\[Rho], 2, 1][t] + 
       1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 3][t] - 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 1][t])},
 {Derivative[1][Subscript[\[Rho], 2, 2]][t] == \[Gamma]/
    2 - \[Gamma] Subscript[\[Rho], 2, 2][t] - 
    I (1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] - 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t]) + 
    1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t]},
 {Derivative[1][Subscript[\[Rho], 2, 3]][t] == 
   1/2 (-\[Gamma] Subscript[\[Rho], 2, 3][
         t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 2, 3][
         t]) - I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 1][t] + 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][
         t] + \[CapitalDelta]b Subscript[\[Rho], 2, 3][t] - 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][t])},
 {Derivative[1][Subscript[\[Rho], 3, 1]][t] == 
   1/2 (-\[Gamma] Subscript[\[Rho], 3, 1][
         t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 1][
         t]) - I (-(1/2) \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][
         t] - 1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 1][
         t] - \[CapitalDelta]a Subscript[\[Rho], 3, 1][t] + 
       1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t])},
 {Derivative[1][Subscript[\[Rho], 3, 2]][t] == 
   1/2 (-\[Gamma] Subscript[\[Rho], 3, 2][
         t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 2][
         t]) - I (-(1/2) \[CapitalOmega]Ra Subscript[\[Rho], 1, 2][
         t] - 1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][
         t] - \[CapitalDelta]b Subscript[\[Rho], 3, 2][t] + 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][t])},
 {Derivative[1][Subscript[\[Rho], 3, 3]][
    t] == -I (-(1/2) \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] - 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] + 
       1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t] + 
       1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][
         t]) - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 3][t]}
}

while the boundary conditions are

{
 {Subscript[\[Rho], 1, 1][0] == 1},
 {Subscript[\[Rho], 1, 2][0] == 0},
 {Subscript[\[Rho], 1, 3][0] == 0},
 {Subscript[\[Rho], 2, 1][0] == 0},
 {Subscript[\[Rho], 2, 2][0] == 0},
 {Subscript[\[Rho], 2, 3][0] == 0},
 {Subscript[\[Rho], 3, 1][0] == 0},
 {Subscript[\[Rho], 3, 2][0] == 0},
 {Subscript[\[Rho], 3, 3][0] == 0}
}

using the numerical values of

{
    \[Gamma] = 0.1;
    
    \[CapitalGamma] = 1;
    
    \[CapitalOmega]Ra = 1;
    
    \[CapitalOmega]Rb = 8;
 
    \[CapitalDelta]a = 1;
    \[CapitalDelta]b = 0.5;
}

These equations are the elements of the Von Neuman equation for three level atomic system. I want to solve the density matrix elements time dependently to see the time dependent behavior of electromagnetically induced transparency.

Improve this question
$\endgroup$
2
  • $\begingroup$ Your ode is complex and the boundary conditions are real valued. Are you looking for real solution? $\endgroup$
    – Ulrich Neumann
    Commented Aug 26, 2022 at 8:33
  • $\begingroup$ @Ulrich Neumann Thanks for your response. I am looking for complex solution. Later on i will use the real and imaginary parts to plot the susceptibility of the system. Regards $\endgroup$
    – Assad Hafiz
    Commented Aug 26, 2022 at 9:15

1 Answer 1

Reset to default
4
$\begingroup$

Edit: Added code to generate 3D plots of selected solution as function of the parameters below

First initialize the constants then flatten the array of equations and initial conditions:

\[Gamma] = 0.1;
\[CapitalGamma] = 1;
\[CapitalOmega]Ra = 1;
\[CapitalOmega]Rb = 8;
\[CapitalDelta]a = 1;
\[CapitalDelta]b = 0.5;

theEqns = {{Derivative[1][Subscript[\[Rho], 1, 1]][
      t] == \[Gamma] (-Subscript[\[Rho], 1, 1][t]) + 
      1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] - 
      I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t]) + \[Gamma]/
      2}, {Derivative[1][Subscript[\[Rho], 1, 2]][
      t] == -(\[Gamma] Subscript[\[Rho], 1, 2][t]) - 
      I (\[CapitalDelta]a Subscript[\[Rho], 1, 2][
           t] - \[CapitalDelta]b Subscript[\[Rho], 1, 2][t] - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 2][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 3][
           t])}, {Derivative[1][Subscript[\[Rho], 1, 3]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 1, 3][
             t]) - \[Gamma] Subscript[\[Rho], 1, 3][t]) - 
      I (\[CapitalDelta]a Subscript[\[Rho], 1, 3][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 2][
           t])}, {Derivative[1][Subscript[\[Rho], 2, 1]][
      t] == -(\[Gamma] Subscript[\[Rho], 2, 1][t]) - 
      I (-(\[CapitalDelta]a Subscript[\[Rho], 2, 1][
             t]) + \[CapitalDelta]b Subscript[\[Rho], 2, 1][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 3][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 1][
           t])}, {Derivative[1][Subscript[\[Rho], 2, 2]][
      t] == \[Gamma] (-Subscript[\[Rho], 2, 2][t]) + 
      1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] - 
      I (1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t]) + \[Gamma]/
      2}, {Derivative[1][Subscript[\[Rho], 2, 3]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 2, 3][
             t]) - \[Gamma] Subscript[\[Rho], 2, 3][t]) - 
      I (\[CapitalDelta]b Subscript[\[Rho], 2, 3][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 1][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] - 
         1/
          2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
           t])}, {Derivative[1][Subscript[\[Rho], 3, 1]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 1][
             t]) - \[Gamma] Subscript[\[Rho], 3, 1][t]) - 
      I (-(\[CapitalDelta]a Subscript[\[Rho], 3, 1][t]) - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 1][
           t])}, {Derivative[1][Subscript[\[Rho], 3, 2]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 2][
             t]) - \[Gamma] Subscript[\[Rho], 3, 2][t]) - 
      I (-(\[CapitalDelta]b Subscript[\[Rho], 3, 2][t]) - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 2][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
           t])}, {Derivative[1][Subscript[\[Rho], 3, 3]][
      t] == -((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 3][
          t]) - I (-(1/
           2) (\[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t]) + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t])}} // Flatten

theInit = {{Subscript[\[Rho], 1, 1][0] == 
     1}, {Subscript[\[Rho], 1, 2][0] == 
     0}, {Subscript[\[Rho], 1, 3][0] == 
     0}, {Subscript[\[Rho], 2, 1][0] == 
     0}, {Subscript[\[Rho], 2, 2][0] == 
     0}, {Subscript[\[Rho], 2, 3][0] == 
     0}, {Subscript[\[Rho], 3, 1][0] == 
     0}, {Subscript[\[Rho], 3, 2][0] == 
     0}, {Subscript[\[Rho], 3, 3][0] == 0}} // Flatten

Then use NDSolveValue to solve the system. Variable solArray contains the nine numerical solutions. First element in solArray, $\texttt{solArray[[1]][t]}$ is the first function $\rho_{1,1}$ and so forth. I use {t,0,1} as an example plotting the real components of the solutions in one plot then a single plot : For example, if you want the value of say $\rho_{2,1}(1/2)$ then use solArray[4][1/2].

solArray = 
  NDSolveValue[
   Join[theEqns, theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho],
     1, 2], Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1], 
    Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3], 
    Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2], 
    Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
Length@solArray
Plot[Re[#[t]] & /@ solArray, {t, 0, 1}, PlotRange -> All]
Plot[Re[solArray[[1]][t]], {t, 0, 1}]

enter image description here

Also, if the array of functions is a bit awkward to deal with, you can individually specify each function as the following and then use them separately (have to use different names on the left though):

{p11, p12, p13, p21, p22, p23, p31, p32, p33} = 
  NDSolveValue[
   Join[theEqns, theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho],
     1, 2], Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1], 
    Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3], 
    Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2], 
    Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
Length@solArray
Plot[Re@p11[t], {t, 0, 1}]

Edit: Added code to write DE as a function of parameters

First write the system as a function of parameters. Now the 9 equations are functions of the parameters of $\texttt{theEqnsF}$:

ClearAll["Global`*"]
theEqns2[\[Gamma]_, \[CapitalGamma]_, \[CapitalOmega]Ra_, \
\[CapitalOmega]Rb_, \[CapitalDelta]a_, \[CapitalDelta]b_] := \
{{Derivative[1][Subscript[\[Rho], 1, 1]][
      t] == \[Gamma] (-Subscript[\[Rho], 1, 1][t]) + 
      1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] - 
      I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t]) + \[Gamma]/
       2}, {Derivative[1][Subscript[\[Rho], 1, 2]][
      t] == -(\[Gamma] Subscript[\[Rho], 1, 2][t]) - 
      I (\[CapitalDelta]a Subscript[\[Rho], 1, 2][
           t] - \[CapitalDelta]b Subscript[\[Rho], 1, 2][t] - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 2][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 3][
           t])}, {Derivative[1][Subscript[\[Rho], 1, 3]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 1, 3][
             t]) - \[Gamma] Subscript[\[Rho], 1, 3][t]) - 
      I (\[CapitalDelta]a Subscript[\[Rho], 1, 3][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 2][
           t])}, {Derivative[1][Subscript[\[Rho], 2, 1]][
      t] == -(\[Gamma] Subscript[\[Rho], 2, 1][t]) - 
      I (-(\[CapitalDelta]a Subscript[\[Rho], 2, 1][
             t]) + \[CapitalDelta]b Subscript[\[Rho], 2, 1][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 3][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 1][
           t])}, {Derivative[1][Subscript[\[Rho], 2, 2]][
      t] == \[Gamma] (-Subscript[\[Rho], 2, 2][t]) + 
      1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] - 
      I (1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t]) + \[Gamma]/
       2}, {Derivative[1][Subscript[\[Rho], 2, 3]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 2, 3][
             t]) - \[Gamma] Subscript[\[Rho], 2, 3][t]) - 
      I (\[CapitalDelta]b Subscript[\[Rho], 2, 3][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 1][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
           t])}, {Derivative[1][Subscript[\[Rho], 3, 1]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 1][
             t]) - \[Gamma] Subscript[\[Rho], 3, 1][t]) - 
      I (-(\[CapitalDelta]a Subscript[\[Rho], 3, 1][t]) - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] + 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 1][
           t])}, {Derivative[1][Subscript[\[Rho], 3, 2]][t] == 
     1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 2][
             t]) - \[Gamma] Subscript[\[Rho], 3, 2][t]) - 
      I (-(\[CapitalDelta]b Subscript[\[Rho], 3, 2][t]) - 
         1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 2][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
           t])}, {Derivative[1][Subscript[\[Rho], 3, 3]][
      t] == -((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 3][
          t]) - I (-(1/2) (\[CapitalOmega]Ra Subscript[\[Rho], 1, 3][
             t]) + 1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t] - 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] + 
         1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t])}} // Flatten

theInit = {{Subscript[\[Rho], 1, 1][0] == 
     1}, {Subscript[\[Rho], 1, 2][0] == 
     0}, {Subscript[\[Rho], 1, 3][0] == 
     0}, {Subscript[\[Rho], 2, 1][0] == 
     0}, {Subscript[\[Rho], 2, 2][0] == 
     0}, {Subscript[\[Rho], 2, 3][0] == 
     0}, {Subscript[\[Rho], 3, 1][0] == 
     0}, {Subscript[\[Rho], 3, 2][0] == 
     0}, {Subscript[\[Rho], 3, 3][0] == 0}} // Flatten

Can now supply the system function $\texttt{theEqnsF}$ to NDSolveValue where now must supply values for the parameters:

   (* theEqns[\[Gamma],\[CapitalGamma],\[CapitalOmega]Ra,\[CapitalOmega]\
Rb,\[CapitalDelta]a,\[CapitalDelta]b] *)
{p11, p12, p13, p21, p22, p23, p31, p32, p33} = 
  NDSolveValue[
   Join[theEqnsF[0.1, 1, 1, 8, 1, 0.5], 
    theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho], 1, 2], 
    Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1], 
    Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3], 
    Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2], 
    Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
Plot[Re@p33[t], {t, 0, 1}]

enter image description here

Now generate a series of plots for $\texttt{p33}$ when $\Delta b$ varies from 0.1 to 5:

plotTable = Table[
   (* theEqns[\[Gamma],\[CapitalGamma],\[CapitalOmega]Ra,\
\[CapitalOmega]Rb,\[CapitalDelta]a,\[CapitalDelta]b] *)
   {p11, p12, p13, p21, p22, p23, p31, p32, p33} = 
    NDSolveValue[
     Join[theEqns[0.1, 1, 1, 8, 1, currentB], 
      theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho], 1, 2], 
      Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1], 
      Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3], 
      Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2], 
      Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
   ParametricPlot3D[{t, currentB, Re@p33[t]}, {t, 0, 1}],
   {currentB, 0.1, 5, 0.1}
   ];
selectedFunction = "p33";
Show[plotTable, BoxRatios -> {1, 1, 1}, PlotRange -> All, 
 AxesLabel -> {Style["t", 14], , Style["\[CapitalDelta]b", 14], 
   Style[selectedFunction, 14]}]

enter image description here

Finally generate a table of values for $\texttt{p33}$ in the form needed for the function ListPlot3D and generate a smooth function of $\texttt{p33}$ over the indicated range:

plotTable = Table[
   (* theEqns[\[Gamma],\[CapitalGamma],\[CapitalOmega]Ra,\
\[CapitalOmega]Rb,\[CapitalDelta]a,\[CapitalDelta]b] *)
   {p11, p12, p13, p21, p22, p23, p31, p32, p33} = 
    NDSolveValue[
     Join[theEqns[0.1, 1, 1, 8, 1, currentB], 
      theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho], 1, 2], 
      Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1], 
      Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3], 
      Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2], 
      Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
           if = p33;
   Table[{t, currentB, if[t]} // N, {t, 0, 1, 1/50}],
   {currentB, 0.1, 5, 0.1}
   ];

ListPlot3D[Flatten[plotTable, 1], BoxRatios -> {1, 1, 1}, 
 AxesLabel -> {Style["t", 14], , Style["\[CapitalDelta]b", 14], 
   Style[selectedFunction, 14]}]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you for your reply and help. can we use different value of [CapitalDelta]b to plot the solution against time and [CapitalDelta]b $\endgroup$
    – Assad Hafiz
    Commented Aug 26, 2022 at 17:46
  • $\begingroup$ @Assad Hafiz: added some code to generate smooth funtion of a selected solution as a function of the parameters. You will need to review the code carefully and adjust it to your needs. $\endgroup$
    – josh
    Commented Aug 26, 2022 at 21:21
  • $\begingroup$ Again thank you for your reply and help. Can we use finite difference method or runge kutta method to solve these equations. $\endgroup$
    – Assad Hafiz
    Commented Aug 27, 2022 at 5:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.