# Visualizing the Schrödinger equation for dipole based interaction Hamiltonian

I want to evaluate and plot the solutions of coupled linear differential equations. For each time frame, I want to display the probability values (corresponding to mod square of these solutions), and have k such frames, going down. Equation is something like this:

sln =
DSolve[
{I c1'[t] == c2[t] + .125 c3[t],
I c2'[t] == c1[t] + c3[t],
I c3'[t] == .125 c1[t] + c2[t],
c1[0] == 1, c2[0] == 0, c3[0] == 0},
{c1[t], c2[t], c3[t]}, t]


When I evaluate the function, it doesn't give me numerical values, as a result I cannot plot the resultant solutions.

Array[Evaluate[{Abs[c1[t]]^2, Abs[c2[t]]^2, Abs[c3[t]]^2}], {t, 0, 20}] /. sln


ArrayPlot doesn't work because the Array doesn't give definite numerical values.

Since,

$\displaystyle{ \sum_{i=1}^{3} |c_i[t]|^2 = 1}$

at all times, a color gradient can be used to represent these values. DensityPlot doesn't help here because it requires two parameters.

• Is this what you are after? ArrayPlot @ Table[{Abs[c1[t]]^2, Abs[c2[t]]^2, Abs[c3[t]]^2} /. sln[[1]], {t, 0, 20}]
– Kuba
Commented Feb 16, 2017 at 11:07
• Thanks, that helps. I can add the RGBColor gradient now. Commented Feb 16, 2017 at 11:10
• @ferro11001 Can you share the source of these equations?
– zhk
Commented Feb 18, 2017 at 13:36

I think, the answer to your question is a straight forward one.

You wanted to plot the solutions of your system, which can be done like this,

{c1sol, c2sol, c3sol} = DSolveValue[{I c1'[t] == c2[t] + .125 c3[t],
I c2'[t] == c1[t] + c3[t], I c3'[t] == .125 c1[t] + c2[t],
c1[0] == 1, c2[0] == 0, c3[0] == 0}, {c1[t], c2[t], c3[t]}, t];

Plot[Evaluate[{Abs[c1sol]^2, Abs[c2sol]^2, Abs[c3sol]^2}], {t, 0, 20},
Frame -> True, PlotStyle -> {Red, Green, Blue}]


Now to extract data from the solution,

Abs[c1sol] /. t -> 0


1

If you want a data table of the solutions, then

TableForm[Table[{t, Abs[c1sol]^2, Abs[c2sol]^2, Abs[c3sol]^2}, {t, 0, 20}],
TableHeadings -> {None, Prepend[{"c1", "c2", "c3"}, t]}]


ArrayPlot[Table[{Abs[c1sol]^2, Abs[c2sol]^2, Abs[c3sol]^2}, {t, 0, 20}],
ColorFunction -> Hue, ColorFunctionScaling -> True]


ListLinePlot[Table[{t, Abs[#]^2}, {t, 0, 20}], Filling -> Axis,
ColorFunction -> Function[{x, y}, RGBColor[1 - y, 0, y]], Frame -> True, Mesh -> Full,
MeshStyle -> {Black, PointSize[Large]}, ImageSize -> 300] & /@ {c1sol, c2sol, c3sol}


{c1sol, c2sol, c3sol} = DSolveValue[{I c1'[t] == c2[t] + .125 c3[t],
I c2'[t] == c1[t] + c3[t], I c3'[t] == .125 c1[t] + c2[t],
c1[0] == c10, c2[0] == c10, c3[0] == c10}, {c1[t], c2[t], c3[t]},
t];
ParametricPlot3D[Evaluate[Table[{Abs[c1sol]^2, Abs[c2sol]^2,
Abs[c3sol]^2} /. {c10 -> m}, {m, 0.0, 2, 0.2}]], {t, 0, 20},
PlotRange -> All, PlotStyle -> {Red, Green, Blue}, PlotPoints -> 100]


• Yes, I got the first graph. I was stuck at writing the simple (now) Table command and getting an ArrayPlot. Physically, this is the Schrödinger equation for dipole based interaction Hamiltonian. Thanks for improving the question as well. Commented Feb 17, 2017 at 2:38
• @NisargBhatt My pleasure. If there is anything else then ask freely?
– zhk
Commented Feb 17, 2017 at 2:53

Maybe something like this !

nmax = 20
psi[n_, t_] = Sin[n t 2 Pi/nmax]^2;
Table[
ListLinePlot[Table[{n, psi[n, t]}, {n, 0, nmax}], Filling -> Axis,
ColorFunction -> Function[{x, y}, RGBColor[1 - y, 0, y]],
Frame -> True, Mesh -> Full, MeshStyle -> {Black, PointSize[Large]},
ImageSize -> 300, PlotLabel -> t]
, {t, 3}]


Or in a compact way

data = Flatten[Table[{n, t, psi[n, t]}, {n, 0, nmax}, {t, 1, 3, 0.1}],1];
ListPlot3D[data, ColorFunction -> Function[{x, y, z}, RGBColor[1 - z, 0, z]],
AxesLabel -> {"n", "t", "psi"}]