Finite difference method for solving coupled differential equations

I am trying to solve the following coupled differential equations: $$\dot{c}_{1}(t)=-\int_0^t \alpha_{1}(t,s)c_{1}(s)+\beta_{1}(t,s)c_{2}(s)\,ds-\frac{i}{\hbar}\lambda c_2(t)e^{+i(\omega_1-\omega_2) t}$$ $$\dot{c}_{2}(t)=-\int_0^t \alpha_{2}(t,s)c_{2}(s)+\beta_{2}(t,s)c_{1}(s)\,ds-\frac{i}{\hbar}\lambda c_1(t)e^{-i(\omega_1-\omega_2) t}$$

I developed a code that uses the finite difference method to solve the equations $$\dot{c}_{1}(t)=-\int_0^t \alpha_{1}(t,s)c_{1}(s)+\beta_{1}(t,s)c_{2}(s)\,ds$$
$$\dot{c}_{2}(t)=-\int_0^t \alpha_{2}(t,s)c_{2}(s)+\beta_{2}(t,s)c_{1}(s)\,ds$$,
where $$\alpha_1(t,s)=\frac{\sqrt{2}\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{01}^2e^{i\omega_{1}(t-s)}}{\hbar^2\left(r_{01}^2+\frac{i\hbar(t-s)}{m_b}\right)^{\frac{3}{2}}}$$,
$$\alpha_2(t,s)=\frac{\sqrt{2}\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{02}^2e^{i\omega_{2}(t-s)}}{\hbar^2\left(r_{02}^2+\frac{i\hbar(t-s)}{m_b}\right)^{\frac{3}{2}}}$$,
$$\beta_1(t,s)=\frac{4\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{01}r_{02}e^{i\omega_{1}(t-s)}m_b^{\frac{3}{2}}e^{i(\omega_1-\omega_2)s}e^{-\frac{\Delta r_c^2m_b}{m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar}}}{{\hbar^2\left(m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar\right)^{\frac{5}{2}}}}$$
and $$\beta_2(t,s)=\frac{4\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{01}r_{02}e^{i\omega_{2}(t-s)}m_b^{\frac{3}{2}}e^{-i(\omega_1-\omega_2)s}e^{-\frac{\Delta r_c^2m_b}{m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar}}}{{\hbar^2\left(m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar\right)^{\frac{5}{2}}}}$$
but I am unsure how to extend this to include the $$-\frac{i}{\hbar}\lambda c_{12}(t)e^{+-i(\omega_1-\omega_2) t}$$ factors (lambda is just a constant). My code to solve the first part of the equations is as follows.

Clear[gab, \[Rho]0, \[HBar], mb, r01, r02, drc, \[Omega]1, \[Omega]2, \
\[Omega]]
drc = 2*r01;
\[HBar] = 1;
\[Omega]1 = 50;
\[Omega]2 = 55;
\[Omega] = (\[Omega]1*\[Omega]2)/(\[Omega]1 + \[Omega]2);
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
r01 = Sqrt[(\[HBar])/(ma*\[Omega]1)];
r02 = Sqrt[(\[HBar])/(ma*\[Omega]2)];
\[Alpha]1[\[Tau]_] := (
Sqrt[2] E^(I (\[Tau]) \[Omega]1) gab^2 \[Pi]^(3/2)
r01^2 \[Rho]0)/(\[HBar]^2 (r01^2 + (I (\[Tau]) \[HBar])/mb)^(
3/2)) ;
\[Alpha]2[\[Tau]_] := (
Sqrt[2] E^(I (\[Tau]) \[Omega]2) gab^2 \[Pi]^(3/2)
r02^2 \[Rho]0)/(\[HBar]^2 (r02^2 + (I (\[Tau]) \[HBar])/mb)^(
3/2)) ;
\[Beta]1[\[Tau]_] := (4 E^(
I  \[Omega]1 \[Tau] - (drc^2 mb)/(
mb (r01^2 + r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(
3/2) \[Pi]^(3/2)
r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar])^(5/2))
\[Beta]2[\[Tau]_] := (4 E^(
I \[Omega]2 \[Tau] - (drc^2 mb)/(
mb (r01^2 + r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(
3/2) \[Pi]^(3/2)
r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar])^(5/2))

dt = 0.015;
nsubint = 1000;
ds = dt/nsubint;
Clear[c1, c2];
c1[0] = 1;
c2[0] = 0;
dp = 30;
\[Lambda] = 10;

Clear[cTtab1, cTtab2];
Do[
corrSum1[n] =
Sum[c1[nn - 1]*
Sum[\[Alpha]1[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[E^(I (\[Omega]1 - \[Omega]2) (m ds)) \[Beta]1[n dt - m ds]*
ds, {m, nsubint (nn - 1), nsubint nn , 1}], {nn, 1, n}];
corrSum2[n] =
Sum[c1[nn - 1]*
Sum[E^(I (\[Omega]2 - \[Omega]1) (m ds)) \[Beta]2[n dt - m ds]*
ds, {m, nsubint (nn - 1), nsubint nn , 1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[\[Alpha]2[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
1}], {nn, 1, n}];
c1[n] = c1[n - 1] - dt*corrSum1[n]-dt*I*\[Lambda]*c2[n - 1];
c2[n] = c2[n - 1] - dt*corrSum2[n]-dt*I*\[Lambda]*c1[n - 1], {n, 1, dp}]
cTtab1 = Table[{n*dt*\[Omega], Abs[c1[n]]^2}, {n, 0, dp}];
cTtab2 = Table[{n*dt*\[Omega], Abs[c2[n]]^2}, {n, 0, dp}];

• It looks like we discussed this problem once on mathematica.stackexchange.com/questions/255697/… Nov 22, 2021 at 4:09
• Hi, thank you very much for your comment! Indeed, we discussed a similar problem, but the introduction of a new function which is only dependent on t makes things more complicated. I thought that creating a new question would be more useful for new people to see it and add there ideas on how to solve it! Nov 22, 2021 at 8:59
• Did you use code from my answer or this is a new code? Nov 22, 2021 at 16:13
• This code is similar to what I posted previously. I finally didn't use the iterative solution you proposed. I managed to make a small adjustment to what I posted to make this working and now I have to include an extra factor, but as this is outside the integral, I am not sure how to. Nov 22, 2021 at 16:20
• Did you test your code or you use it as it is? In my answer there are 3 methods are compared to test your code. Nov 23, 2021 at 3:36

We can solve this problem using FDM and implicit method of integration differential equations as follows:

Clear[gab, \[Rho]0, \[HBar], mb, r01, r02, drc, \[Omega]1, \[Omega]2, \
\[Omega]]
drc = 2*r01;
\[HBar] = 1; \[Lambda] = 10; k = I \[Lambda]/\[HBar];
\[Omega]1 = 50;
\[Omega]2 = 55; dom = \[Omega]2 - \[Omega]1;
\[Omega] = (\[Omega]1*\[Omega]2)/(\[Omega]1 + \[Omega]2);
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
r01 = Sqrt[(\[HBar])/(ma*\[Omega]1)];
r02 = Sqrt[(\[HBar])/(ma*\[Omega]2)];
f1[\[Tau]_] := (Sqrt[
2] E^(I (\[Tau]) \[Omega]1) gab^2 \[Pi]^(3/
2) r01^2 \[Rho]0)/(\[HBar]^2 (r01^2 + (I (\[Tau]) \[HBar])/
mb)^(3/2));
f2[\[Tau]_] := (Sqrt[
2] E^(I (\[Tau]) \[Omega]2) gab^2 \[Pi]^(3/
2) r02^2 \[Rho]0)/(\[HBar]^2 (r02^2 + (I (\[Tau]) \[HBar])/
mb)^(3/2));
g1[\[Tau]_] := (4 E^(I \[Omega]1 \[Tau] - (drc^2 mb)/(mb (r01^2 +
r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
2) r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar])^(5/2));
g2[\[Tau]_] := (4 E^(I \[Omega]2 \[Tau] - (drc^2 mb)/(mb (r01^2 +
r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
2) r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) +
2 I (\[Tau]) \[HBar])^(5/2));

dt = 0.01;
nsubint = 10;
ds = dt/nsubint;
Clear[c1, c2];
c1[0] = 1; c1[-1] = 1;
c2[0] = 0; c2[-1] = 0;

Clear[cTtab1, cTtab2];
Do[corrSum1[n] =
Sum[c1[nn - 1]*
Sum[f1[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[g1[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}];
corrSum2[n] =
Sum[c1[nn - 1]*
Sum[g2[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[f2[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}];
c1[n] = -((-c1[-1 + n] + 1/2 dt E^(I dom dt (-1 + n)) k c2[-1 + n] +
dt corrSum1[n] -
1/2 dt E^(I dom dt n)
k (1/2 dt E^(-I dom dt (-1 + n)) k c1[-1 + n] - c2[-1 + n] +
dt corrSum2[n]))/(1 - (dt^2 k^2)/4));
c2[n] = -(1/(-4 + dt^2 k^2)) E^(-I dom dt (-1 + n) -
I dom dt n) (-2 dt E^(I dom dt (-1 + n)) k c1[-1 + n] -
2 dt E^(I dom dt n) k c1[-1 + n] +
4 E^(I dom dt (-1 + n) + I dom dt n) c2[-1 + n] +
dt^2 E^(2 I dom dt (-1 + n)) k^2 c2[-1 + n] +
2 dt^2 E^(I dom dt (-1 + n)) k corrSum1[n] -
4 dt E^(I dom dt (-1 + n) + I dom dt n) corrSum2[n]);, {n, 0,
500}]


Visualization

cTtab1 = Table[{(n + 1)*dt, Abs[c1[n]]}, {n, 0, 500}];
cTtab2 = Table[{(n + 1)*dt, Abs[c2[n]]}, {n, 0, 500}];

pl1 = ListLinePlot[{cTtab1, cTtab2}, PlotRange -> All,
PlotLegends -> {"c1", "c2"}, AxesLabel -> {"t", ""}]


We can compare FDM result (points) with iterative method (solid lines) explained in my answer here

X[0][t_] := 1;
Y[0][t_] := 0; ds1 =
1/100; tmax = 1; nmax = 10; Do[{X[i], Y[i]} =
NDSolveValue[{x'[t]/
tmax == -tmax t ds1/2 Sum[
f1[tmax (t - t/2 (s + 1))] X[i - 1][t/2 (s + 1)] +
g1[tmax (t - t/2 (s + 1))] Y[i - 1][t/2 (s + 1)], {s, -1, 1,
ds1}] - k y[t] Exp[I dom tmax t],
y'[t]/tmax == -tmax t ds1/2 Sum[
g2[tmax (t - t/2 (s + 1))] X[i - 1][t/2 (s + 1)] +
f2[tmax (t - t/2 (s + 1))] Y[i - 1][t/2 (s + 1)], {s, -1, 1,
ds1}] - k x[t] Exp[-I dom tmax t], x[0] == 1,
y[0] == 0}, {x, y}, {t, 0, 1}];, {i, 1, nmax}]

pl2 = Plot[Evaluate[{X[nmax][t], Y[nmax][t]} // Abs], {t, 0, 1},
PlotLegends -> {"c1", "c2"}];
pl3 = ListPlot[{Take[cTtab1, 100], Take[cTtab2, 100]},
PlotRange -> All];
Show[pl2, pl3]


• Hi, thank you for the great answer. I have some doubts. In the FDM, how have you accounted for the exp(i dom s) factor? And what is the motivation behind using c1[n]=-... and not just the positive? Also, is there a way to extend the iterative method beyond t=1.0? Nov 28, 2021 at 12:28
• To get c1[n],c2[n] just evaluate Solve[c1[n] == c1[n - 1] - dt*corrSum1[n] - (c2[n - 1] Exp[I dom (n - 1) dt] + c2[n] Exp[I dom n dt])/2 dt k && c2[n] == c2[n - 1] - dt*corrSum2[n] - (c1[n - 1] Exp[-I dom (n - 1) dt] + c1[n] Exp[-I dom n dt])/2 dt k, {c1[n], c2[n]}]. Nov 28, 2021 at 12:39
• To extend iterative solution up to tmax=3 we can put nmax=25 and use FDM solution as starting point. Nov 28, 2021 at 17:44
• Also we can use as starting point X[0][t_] := Exp[-t/tmax]; Y[0][t_] := 0; Nov 28, 2021 at 17:54