# Crank-Nicolson scheme for the Schrödinger equation

I am trying to solve the 1-D Schrödinger equation with the Crank-Nicholson method. I use the atomic units:

• time < 10^-14 => 413 a.u.t.
• length 30 * 10^-9 m => 567 Bohr radii

a = 50;
σ = a/10;
p[x_, 0] := Exp[-(x^2)/(2 σ^2)]; (* граничные условия *)
p[\[PlusMinus]a, t_] := 0; (* н.у. *)
NN = 100;
h = 1/NN;
repl = x -> ((i*6*a)/NN - 3*a);

eq = D[Subscript[p, i][t], t] == I/2(Subscript[p, i - 1][t] + 2 Subscript[p, i][t] -
Subscript[p, i + 1][t])/h^2(*+D[u[x],x]*Subscript[p,i][t]*)/. repl;

Table[eq, {i, 0, NN}];
boundary = {Subscript[p, 0][t] == 0, Subscript[p, NN][t] == 0};
p0[x_] := Exp[-(x^2)/(2 σ^2)];
Cauchy = Table[Subscript[p, i][0] == p0[x] /. repl, {i, 1, NN - 1}]; (*задача Коши*)

eqns = Join[Table[eq, {i, 1, NN - 1}], boundary, Cauchy];
sol = NDSolve[N[eqns],
N[Table[Subscript[p, i][t], {i, 0, NN}]], {t, 0, 413}];


and get this error:

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. >>

I want to get something like this:

sol = NDSolve[{I D[u[t, x], t] == (-1/2)D[u[t, x], {x, 2}],
u[0., x] == Exp[-(x^2./(2*σ^2))], u[t, a] == 0,
u[t, -a] == 0}, u, {t, 0, 4130}, {x, -a, a}, (*MaxStepSize->0.01,*)
AccuracyGoal -> 3, PrecisionGoal -> 3];

Animate[Plot[Evaluate[Abs[u[t, x] /. First[sol]]^2], {x, -a, a},
PlotRange -> {0, 1}], {t, 0, 413}]


What's wrong with my code and how to fix it?

-
Back in the late 90s I recall Ron Knapp writing an example for CN for Schrodinger. Please search the Wolfram library and I think you should be able to locate it. – Mike Honeychurch Dec 13 '13 at 2:18
@MikeHoneychurch I tried to search but can't find Ron Knapp's example, could you search if you have time? I'm very interested in this. But I find this example by Terry Robb, which solves the TDSE using a external Fortran subroutine SCHROED. And SCHROED uses a Crank-Nicholson scheme. – xslittlegrass Mar 6 '14 at 4:51