# Crank-Nicolson scheme for Schrodinger equation

I try to solve 1D Schrodinger equation with Crank-Nickolson methode. I use atomic units: time < 10^-14 => 413 a.u.t. length 30 * 10^-9 m => 567 Bohr radii

a = 50;
\[Sigma] = a/10;
p[x_, 0] := Exp[-(x^2)/(2 \[Sigma]^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 \[Sigma]^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 have 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*\[Sigma]^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 at 4:51