Solving coupled ODEs with NDSolve

Question

Given the system of differential equations:

$i\frac{d}{dt}c_{1}(t)=H_{12} E(t)e^{(iw_{12}t)}c_2(t)+H_{13} E(t)e^{(iw_{13}t)}c_3(t)$, $i\frac{d}{dt}c_{2}(t)=H_{21} E(t)e^{(iw_{21}t)}c_2(t)+H_{13} E(t)e^{(iw_{13}t)}c_3(t)$, $i\frac{d}{dt}c_{3}(t)=H_{31} E(t)e^{(iw_{31}t)}c_2(t)+H_{32} E(t)e^{(iw_{32}t)}c_2(t)$,

where $H_{kn}=0.005, w_{12}=-0.001, w_{21}=0.001,w_{13}=1.000,w_{31}=-1.000,w_{32}=0.999,w_{23}=0.999$ $E(t)=E_0 sech(\frac{t-400}{150})^2 sin(t)$, where $E_0=9.0$,

I wish to plot $|c_1(t)|^2,|c_2(t)|^2,|c_3(t)|^2 $vs $t$,

I'm relatively new in using Mathematica. I tried to solve this differential equation by using NDSolve. Although I get the output as an interpolating function, it still does not generate any plot. Any help will be highly helpful!

My code is

Clear[x, y, z, s, c1, c2, c3, t, w12, w21, w23, w23, w32, w13, w31, e0, ham]
w12 = -0.001;
w21 = 0.001;
w13 = -1.000;
w31 = 1.000;
w23 = -0.999;
w32 = 0.999;
ham = 0.005;
e0 = 9.0;
paramND = {
  I*D[c1[t], t] == 
    ham*(e0*Sech[(t - 400)/150]^2*Sin[t])*((Exp[I*w12*t]*c2[t]) + (Exp[I*w13*t]*c3[t])),
  I*D[c2[t], t] == 
    ham*(e0*Sech[(t - 400)/150]^2Sin[t])*((Exp[I*w21*t]*c1[t])+ (Exp[I*w23*t]*c3[t])),
  I*D[c3[t], t] == 
    ham*(e0*Sech[(t - 400)/150]^2*Sin[t])*((Exp[I*w31*t]*c1[t]) + (Exp[I*w32*t]*c2[t])),
   c1[0.01] == 1.0, c2[0.01] == 0.0, c3[0.01] == 0.0};
numSol = 
  NDSolve[paramND, {c1[t], c2[t], c3[t]}, {t, 0.01, 1}, 
    MaxSteps -> ∞, Method -> "ExplicitRungeKutta"] ;
paramPlot = {c1[t]^2, c2[t]^2, c3[t]^2} //. numSol;
Plot[Evaluate @ paramPlot, {t, 0.01, 800},
  PlotRange -> {0.01,1.0}, PlotStyle -> {{Black, Thick}}]

Answer 1

it works, you just have complex functions, also the range you had was not helping. It is also better I think to use {c1, c2, c3} instead of {c1[t], c2[t], c3[t]}

Clear[x, y, z, s, c1, c2, c3, t, w12, w21, w23, w23, w32, w13, w31, e0, ham]
w12 = -0.001;
w21 = 0.001;
w13 = -1.000;
w31 = 1.000;
w23 = -0.999;
w32 = 0.999;
ham = 0.005;
e0 = 9.0;
paramND = {I*D[c1[t], t] == ham*(e0*Sech[(t - 400)/150]^2*Sin[t])*
   ((Exp[I*w12*t]*c2[t]) + (Exp[I*w13*t]*c3[t])), 
   I*D[c2[t], t] == ham*(e0*Sech[(t - 400)/150]^2 Sin[t])*((Exp[I*w21*t]*c1[t]) + 
  (Exp[I*w23*t]*c3[t])), 
   I*D[c3[t], t] == ham*(e0*Sech[(t - 400)/150]^2*Sin[t])*((Exp[I*w31*t]*c1[t]) + 
   (Exp[I*w32*t]*c2[t])), c1[0.01] == 1.0, 
   c2[0.01] == 0.0, c3[0.01] == 0.0};
numSol = First@NDSolve[paramND, {c1, c2, c3}, {t, 0.01, 1}, MaxSteps -> \[Infinity], 
    Method -> "ExplicitRungeKutta"];
paramPlot = c3 /. numSol;
Plot[Re@paramPlot[t], {t, 0.01, 800}, PlotRange -> All, PlotStyle -> {{Black, Thick}}]

Plot[{Re@paramPlot[t], Im@paramPlot[t]}, {t, 0.01, 1}, PlotRange -> All, 
 PlotTheme -> "Detailed"]

Update for comment:

Plot[{Abs[(c1 /. numSol)[t]]^2, Abs[(c2 /. numSol)[t]]^2, 
 Abs[(c3 /. numSol)[t]]^2}, {t, 0.01, 800}, PlotRange -> All, PlotTheme -> "Detailed"]