# Finite difference method to solve coupled differential equations

I am trying to solve the following coupled differential equations: $$\dot{c}_1(t)=-\int_0^tf(t-s)c_1(s)+g(t-s)c_2(s)ds$$, $$\dot{c}_2(t)=-\int_0^tg(t-s)c_1(s)+f(t-s)c_2(s)ds$$, where $$f(t-s)=\frac{\sqrt{2}\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_0^2e^{i\omega(t-s)}}{\hbar^2\left(r_0^2+\frac{i\hbar(t-s)}{m_b}\right)^{\frac{3}{2}}}$$ and $$g(t-s)=\frac{\sqrt{2}e^{i\omega(t-s)-\frac{\Delta r_c^2m_b}{2m_br_0^2+2i(t-s)\hbar}}g_{ab}^2m_b^{\frac{3}{2}}\pi^{\frac{3}{2}}r_0^2\rho_0\left(-\Delta r_c^2m_b+m_br_0^2+i(t-s)\hbar\right)}{\hbar^2\left(m_br_0^2+i\hbar(t-s)\right)^{\frac{5}{2}}}$$.

I developed a code that uses the finite difference method to solve an equation in the form $$\dot{c}_1(t)=-\int_0^tf(t-s)c_1(s)ds$$, but I am unsure how to extend this to solve the coupled differential equations. My attempt to do so is as follows:

Clear[gab, \[Rho]0, \[HBar], \[Omega]0, mb, \[Tau], k, r0, t, s, drc]
drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
\[Omega]0 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (
Sqrt[2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/2)
r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/mb)^(
3/2)) ;
g[\[Tau]_] := (Sqrt[2] E^(
I \[Omega]0 \[Tau] - (drc^2 mb)/(
2 mb r0^2 + 2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/2)
r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 +
I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 +
I (\[Tau]) \[HBar])^(5/2));
dt = 0.1;
nsubint = 1000;
ds = dt/nsubint;
Ttab = Table[T, {T, 0, 10, dt}];
Stab = Table[s, {s, 0, dt - ds, ds}];
(*Ti[T_]:=Position[Ttab,T][[1]]*)
Clear[Ti];
Clear[c1, c2];
c1[0] = 1;
c2[0] = 0;
Clear[cTtab1, cTab2];
(*c1[1]=c1[0]-dt(c1[0]*f[dt]+c2[0]*g[dt];*)
\
(*c2[1]=c2[0]-dt(c1[0]*g[dt]+c2[0]*f[dt];*)
Do[
corrSum1[n] =
Sum[c1[nn - 1]*
Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
1}], {nn, 1, n}];
corrSum2[n] =
Sum[c1[nn - 1]*
Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
1}], {nn, 1, n}];
c1[n] = c1[n - 1] - dt*corrSum1[n];
c2[n] = c2[n - 1] - dt*corrSum2[n], {n, 0, 100}]
cTtab1 = Table[{n*dt, c1[n]}, {n, 0, 100}]
cTtab2 = Table[{n*dt, c2[n]}, {n, 0, 100}]
FDiff1 = ListPlot[Abs[cTtab1], PlotStyle -> Orange,
PlotLegends -> {"Finite Difference c1"}]
FDiff2 = ListPlot[Abs[cTtab2], PlotStyle -> Orange,
PlotLegends -> {"Finite Difference c2"}]


The absolute value of c1 and c2 should go between 0 and 1.

• Can you add the definition of g[tau] to your code? Sep 16, 2021 at 21:00
• You're right, just added it! Sep 17, 2021 at 12:44
• This is a system of integrodifferential equations. Since f, g are complex we have 4 real equations as it shown by Cesario. Sep 18, 2021 at 8:16
• @JoseEnriqueAroca Your code is quite well, but it is not so precise as, for instance, iterative method (see my answer). Sep 20, 2021 at 8:52

We can use the original code after small correction as follows

Clear[gab, \[Rho]0, \[HBar], \[Omega]0, mb, \[Tau], k, r0, t, s, drc]
drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
\[Omega]0 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (Sqrt[
2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/
2) r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/
mb)^(3/2));
g[\[Tau]_] := (Sqrt[
2] E^(I \[Omega]0 \[Tau] - (drc^2 mb)/(2 mb r0^2 +
2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
2) r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 +
I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 +
I (\[Tau]) \[HBar])^(5/2));
dt = 0.01;
nsubint = 10;
ds = dt/nsubint;
Ttab = Table[T, {T, 0, 10, dt}];
Stab = Table[s, {s, 0, dt - ds, ds}];
(*Ti[T_]:=Position[Ttab,T][[1]]*)
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[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}];
corrSum2[n] =
Sum[c1[nn - 1]*
Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}];
c1[n] = c1[n - 1] - dt*corrSum1[n];
c2[n] = c2[n - 1] - dt*corrSum2[n], {n, 0, 100}]

cTtab1 = Table[{n*dt, Abs[c1[n]]}, {n, 0, 100}];
cTtab2 = Table[{n*dt, Abs[c2[n]]}, {n, 0, 100}];

FDiff1 = ListPlot[cTtab1, PlotStyle -> Orange,
PlotLegends -> {"Finite Difference c1"}]
FDiff2 = ListPlot[cTtab2, PlotStyle -> Orange,
PlotLegends -> {"Finite Difference c2"}]


The second method is iteration with substitution s->t/2 (s+1), and therefore Volterra equations are reduced to Fredholm equations as follows

X[0][t_] := 1;
Y[0][t_] := 0; ds =
1/100; nmax = 15; Do[{X[i], Y[i]} =
NDSolveValue[{x'[
t] == -t ds/2 Sum[
f[t - t/2 (s + 1)] X[i - 1][t/2 (s + 1)] +
g[t - t/2 (s + 1)] Y[i - 1][t/2 (s + 1)], {s, -1, 1, ds}],
y'[t] == -t ds/2 Sum[
g[t - t/2 (s + 1)] X[i - 1][t/2 (s + 1)] +
f[t - t/2 (s + 1)] Y[i - 1][t/2 (s + 1)], {s, -1, 1, ds}],
x[0] == 1, y[0] == 0}, {x, y}, {t, 0, 1}];, {i, 1, nmax}]


Visualization of FDM and iterative solution

{Show[Plot[Evaluate[Abs[X[nmax][t]]], {t, 0, 1},
FrameLabel -> {t, "c1"}, Frame -> True], FDiff1],
Show[Plot[Evaluate[Abs[Y[nmax][t]]], {t, 0, 1},
FrameLabel -> {t, "c2"}, Frame -> True], FDiff2]}


We also can compare the iterative method with code @Cesario. After small modification we have

Clear["Global*"]
n = 100;
tmax = 1;
dT = tmax/n; drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
\[Omega]0 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (Sqrt[
2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/
2) r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/
mb)^(3/2));
g[\[Tau]_] := (Sqrt[
2] E^(I \[Omega]0 \[Tau] - (drc^2 mb)/(2 mb r0^2 +
2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
2) r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 +
I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 +
I (\[Tau]) \[HBar])^(5/2));
fr[k_] := Re[N[f[k dT]]]
fi[k_] := Im[N[f[k dT]]]
gr[k_] := Re[N[g[k dT]]]
gi[k_] := Im[N[g[k dT]]]

vars = Flatten[Table[{c1r[k], c1i[k], c2r[k], c2i[k]}, {k, 0, n}]];
equ1r[k_] := (c1r[k] - c1r[k - 1])/dT +
dT Sum[(-c1i[j] fi[1 - j + k] + c1r[j] fr[1 - j + k] -
c2i[j] gi[1 - j + k] + c2r[j] gr[1 - j + k]), {j, 1, k}]
equ1i[k_] := (c1i[k] - c1i[k - 1])/dT +
dT Sum[(c1r[j] fi[1 - j + k] + c1i[j] fr[1 - j + k] +
c2r[j] gi[1 - j + k] + c2i[j] gr[1 - j + k]), {j, 1, k}]
equ2r[k_] := (c2r[k] - c2r[k - 1])/dT +
dT Sum[(-c2i[j] fi[1 - j + k] + c2r[j] fr[1 - j + k] -
c1i[j] gi[1 - j + k] + c1r[j] gr[1 - j + k]), {j, 1, k}]
equ2i[k_] := (c2i[k] - c2i[k - 1])/dT +
dT Sum[(c2r[j] fi[1 - j + k] + c2i[j] fr[1 - j + k] +
c1r[j] gi[1 - j + k] + c1i[j] gr[1 - j + k]), {j, 1, k}]
equs = Flatten[
Table[{equ1r[k] == 0, equ1i[k] == 0, equ2r[k] == 0,
equ2i[k] == 0}, {k, 1, n}]];
equstot =
Join[equs, {c1r[0] - 1 == 0, c1i[0] == 0, c2r[0] == 0, c2i[0] == 0}];
sol = Solve[equstot, vars][[1]];

p1 = ListPlot[Table[{k dT, Norm[{c1r[k], c1i[k]}] /. sol}, {k, 0, n}],
PlotStyle -> Red]

p2 = ListPlot[Table[{k dT, Norm[{c2r[k], c2i[k]}] /. sol}, {k, 0, n}],
PlotStyle -> Red]


Comparison of iterative method and Cesario's method (X, Y are taken from the code above)

{Show[Plot[Evaluate[Abs[X[nmax][t]]], {t, 0, 1},
FrameLabel -> {t, "c1"}, Frame -> True], p1],
Show[Plot[Evaluate[Abs[Y[nmax][t]]], {t, 0, 1},
FrameLabel -> {t, "c2"}, Frame -> True], p2]}


Note, that all 3 methods give us qualitatively same results, but they are differ in detail.
For parametric research we can use Module, for example,

plot[om_] := Module[{\[Omega]0 = om, f, g}, drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (Sqrt[
2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/
2) r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/
mb)^(3/2));
g[\[Tau]_] := (Sqrt[
2] E^(I \[Omega]0 \[Tau] - (drc^2 mb)/(2 mb r0^2 +
2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
2) r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 +
I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 +
I (\[Tau]) \[HBar])^(5/2));
dt = 0.01;
nsubint = 10;
ds = dt/nsubint;
Ttab = Table[T, {T, 0, 10, dt}];
Stab = Table[s, {s, 0, dt - ds, ds}];
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[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[
g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 1}], {nn,
1, n}];
corrSum2[n] =
Sum[c1[nn - 1]*
Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}] +
Sum[c2[nn - 1]*
Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn,
1}], {nn, 1, n}];
c1[n] = c1[n - 1] - dt*corrSum1[n];
c2[n] = c2[n - 1] - dt*corrSum2[n], {n, 0, 100}];
cTtab1 = Table[{n*dt, Abs[c1[n]]}, {n, 0, 100}];
cTtab2 = Table[{n*dt, Abs[c2[n]]}, {n, 0, 100}]; {cTtab1, cTtab2}]


With this code we can plot solution as function of $$\omega 0$$ as follows

pl = Table[
ListPlot[plot[om], PlotRange -> All, PlotLegends -> {"c1", "c2"},
PlotLabel -> Row[{"\[Omega]0 = " , om}],
AxesLabel -> {"t", ""}], {om, {5, 10, 15}}]


We can also combine curves in two plots

pl1 = ListPlot[Table[plot[om][[1]], {om, {5, 10, 15}}],
PlotRange -> All, AxesLabel -> {"t", "c1"},
PlotLegends -> {5, 10, 15}]

pl2 = ListPlot[Table[plot[om][[2]], {om, {5, 10, 15}}],
PlotRange -> All, AxesLabel -> {"t", "c2"},
PlotLegends -> {5, 10, 15}]


• Hi, thank you very much for your detailed answer. I have a small question. Say I want to use the original code to do some plots with changing parameters (e.g. different Omega_1 values) and plot the results on the same axis. What would be the most efficient way to do this? Can I create an array of omega_2 values or something similar? Sep 24, 2021 at 11:15
• @JoseEnriqueAroca Your code not depends on $\omega 1$, but depends on $\omega 0$. We can use your code as Module to plot solution vs different $\omega 0$ - see update to my answer. Sep 24, 2021 at 16:44
• Hi yeah, sorry for the typo, I meant 𝜔0. Is there any way to show the plots for the three 𝜔0 values on the same axis (i.e show the 6 lines on the same axis) Sep 25, 2021 at 15:57
• @JoseEnriqueAroca May be it could better to combine curves in two plots? See update to my answer. Sep 26, 2021 at 3:31

To have a coarse view of $$c_1,c_2$$ we can proceed as follows

Clear["Global*"]
n = 100;
tmax = 5;
dT = tmax/n;
rc = 1;
\[HBar] = 1;
\[Omega] = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
r0 = 1;
\[CapitalDelta] = 1;
f[k_] := (E^(I k dT \[Omega]) gab^2 \[Pi]^(3/2) r0^2\[Rho]0)/(Sqrt[2] \[HBar]^2 (r0^2 + (I k dT \[HBar])/mb)^(3/2));
g[k_] := (Sqrt[2] E^(I k dT \[Omega] - \[CapitalDelta] rc^2 mb/(2 mb r0^2 + 2 I k dT \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/2) r0^2 \[Rho]0) (-\[CapitalDelta] rc^2 mb + mb r0^2 + I k dT \[HBar])/( \[HBar]^2 (mb r0^2 + (I k dT \[HBar]))^(5/2));

fr[k_] := Re[N[f[k]]]
fi[k_] := Im[N[f[k]]]
gr[k_] := Re[N[g[k]]]
gi[k_] := Im[N[g[k]]]

vars = Flatten[Table[{c1r[k], c1i[k], c2r[k], c2i[k]}, {k, 0, n}]];
equ1r[k_] := (c1r[k] - c1r[k - 1])/dT + dT Sum[(- c1i[j] fi[1 - j + k] + c1r[j] fr[1 - j + k] - c2i[j] gi[1 - j + k] + c2r[j] gr[1 - j + k]), {j, 1, k}]
equ1i[k_] := (c1i[k] - c1i[k - 1])/dT + dT Sum[( c1r[j] fi[1 - j + k] + c1i[j] fr[1 - j + k] + c2r[j] gi[1 - j + k] + c2i[j] gr[1 - j + k]), {j, 1, k}]
equ2r[k_] := (c2r[k] - c2r[k - 1])/dT + dT Sum[(- c2i[j] fi[1 - j + k] + c2r[j] fr[1 - j + k] - c1i[j] gi[1 - j + k] + c1r[j] gr[1 - j + k]), {j, 1, k}]
equ2i[k_] := (c2i[k] - c2i[k - 1])/dT + dT Sum[( c2r[j] fi[1 - j + k] + c2i[j] fr[1 - j + k] + c1r[j] gi[1 - j + k] + c1i[j] gr[1 - j + k]), {j, 1, k}]
equs = Flatten[Table[{equ1r[k], equ1i[k], equ2r[k], equ2i[k]}, {k, 1, n}]];
equstot = Join[equs, {c1r[0] - 1/2, c1i[0], c2r[0], c2i[0] - 1/2}];
sol = Solve[equstot == 0, vars][[1]];

(* ListPlot[Table[{k dT, c1r[k]}, {k, 0, n}] /. sol]
ListPlot[Table[{k dT, c1i[k]}, {k, 0, n}] /. sol] *)
ListPlot[Table[{c1r[k], c1i[k]}, {k, 0, n}] /. sol]
(* ListPlot[Table[{k dT, c2r[k]}, {k, 0, n}] /. sol]
ListPlot[Table[{k dT, c2i[k]}, {k, 0, n}] /. sol] *)
ListPlot[Table[{c2r[k], c2i[k]}, {k, 0, n}] /. sol]


• Hi, thank you very much for your answer. I think that this would work fine. However, I have been asked to explicitly use the finite difference method shown above. I have edited the code in the question with my attempt to include c1 and c2 in the same Do Loop, but for some reason, I get answers for c1 and c2 bigger than 1 (they should never exceed this value). You can take a look if you like. Sep 17, 2021 at 12:46
• Nice code (+1), but you have typo in definition of fr[k_] := Re[N[f[k]]]; fi[k_] := Im[N[f[k]]]; gr[k_] := Re[N[g[k]]]; gi[k_] := Im[N[g[k]]]. It should be fr[k_] := Re[N[f[k dT]]];... Please, have a look Figure 3 in my answer with comparison your method and iterative method. Sep 20, 2021 at 8:50
• Thanks for your appreciation (+1). In my case, the k*dT transformation is included into the function definition. Sep 20, 2021 at 9:30