Time dependent Schrödinger equation in 2D

Question

I have the following Schrödinger equation in $2D$:

\begin{cases} \partial_t \Psi(x,t) = V(x,t) \Psi(x,t) \quad x \in [-10,10]^2\\ \Psi(x,0)=\exp( \frac{1}{2} (-x^2+y^2)) \end{cases}

where the potential $V(x,t)=\mathbb{i} \Bigl( \frac{1}{2} \Delta - (x^2+y^2) - \sin^2(t) (x+y) \Bigr)$ with homogeneous Dirichlet boundary conditions. I need the solution at time $T=1$.

Using second order finite differences, I obtain the following plot, plotting $|U|$ at $T=1$:

with the following colormap

I'd like to use Mathematica to check my results and to try what comes out by changing some parameters, but I don't know how to solve it properly. Could someone show the plot of the surface I should obtain with Mathematica, and, if possible, the right code-snippet?

EDIT:

I had a different initial data, now my plot seems to agree with the on of Henrik

Answer 1

Something like the following should do. It employ the finite element method.

Ω =   DiscretizeRegion[Rectangle[{-10, -10}, {10, 10}], MaxCellMeasure -> (1 -> 0.5)];
sol = NDSolveValue[
   {
    D[Ψ[x, y, t], t] == I/2 Laplacian[Ψ[x, y, t], {x, y}] - I ((x^2 + y^2) + (x + y) Sin[t]^2) Ψ[x, y, t], 
    DirichletCondition[Ψ[x, y, t] == 0, True],
    Ψ[x, y, 0] == Exp[-1/2 (x^2 + y^2)]
    },
   Ψ,
   {t, 0, 1},
   {x, y} ∈ Ω
   ];
Plot3D[Abs[sol[x, y, 1]], {x, y} ∈ Ω, PlotRange -> All, AxesLabel -> {"x", "y", "|Ψ|"}]

Looks a bit different from OP's solution, but that could be to a copying error... Anyways, this shows roughly how the PDE can be solved.

For further details (in particular on how to increase the accuracy of the solution), please refer to the documentation (https://reference.wolfram.com/language/FEMDocumentation/tutorial/FiniteElementOverview.html).

Finding the maximum:

NMaximize[{Abs[sol[x, y, 1]], -10 <= x <= 10, -10 <= y <= 10}, {x, y}]

{1.38754, {x -> -0.0632606, y -> -0.0637582}}

Answer 2

FiniteElement isn't necessary for this problem. The old good TensorProductGrid handles the problem quite well:

system = With[{Ψ = Ψ[x, y, t]}, 
          {D[Ψ, t] == I (Laplacian[Ψ, {x, y}]/2 - ((x^2 + y^2) + Sin[t]^2 (x + y)) Ψ),
           Ψ == 0 /. {{x -> -10}, {x -> 10}, {y -> -10}, {y -> 10}},
           Ψ == Exp[-1/2 (x^2 + y^2)] /. t -> 0}];

sol = NDSolveValue[system, Ψ, {t, 0, 1}, {x, -10, 10}, {y, -10, 10}];

Plot3D[Abs@sol[x, y, 1], {x, -10, 10}, {y, -10, 10}, PlotRange -> All, PlotPoints -> 50]

NMaximize[Abs[sol[x, y, 1]], {x, y}]   
(* {1.4014, {x -> -0.0593488, y -> -0.0593488}} *)

Test passes in v12.1.1.

---

Futher tests show v9.0.1 and v8.0.4 have difficulty in solving the system with defaullt setting, so this turns out to be another example indicating NDSolve is improved silently these years. Nevertheless, with the magic of Pseudospectral, we can still solve the problem in v8 and v9:

If[$VersionNumber < 9, Laplacian = D[#, x, x] + D[#, y, y] &;
  NDSolveValue = #2 /. First@NDSolve[##] &];

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

system = With[{Ψ = Ψ[x, y, t]}, 
          {D[Ψ, t] == I (Laplacian[Ψ, {x, y}]/2 - ((x^2 + y^2) + Sin[t]^2 (x + y)) Ψ),
           Ψ == 0 /. {{x -> -10}, {x -> 10}, {y -> -10}, {y -> 10}},
           Ψ == Exp[-1/2 (x^2 + y^2)] /. t -> 0}];

sol = NDSolveValue[system, Ψ, {t, 0, 1}, {x, -10, 10}, {y, -10, 10}, 
    Method -> mol[55]]; // AbsoluteTiming
(* v8.0.4: {178.4673377, Null} *)
(* v9.0.1: {40.305892, Null} *)

FindMaximum[Abs@sol[x, y, 1], {x, y}]
(* v8.0.4: {1.38975, {x -> -0.0438577, y -> -0.0438577}} *)
(* v9.0.1: lstol warning, {1.38918, {x -> -0.0439239, y -> -0.043924}} *)

NMaximize isn't used to find the maximum because it spits out a Experimental`NumericalFunction[…] as output in v8 and v9, which is obviously a (now fixed) bug.

Answer 3

You can simply solve this equation using NDSolve.

Note, I rewrote your equation a bit more towards standard form.

V[x_, y_, t_] := (x^2 + y^2 +  Sin[t]^2 (x + y));
eq = {I  Derivative[0, 0, 1][f][x, y, 
      t] == -Laplacian[f[x, y, t], {x, y}]/2 + V[x, y, t] f[x, y, t], 
   f[x, y, 0] == Exp[-1/2 (x^2 + y^2)], 
   DirichletCondition[f[x, y, t] == 0, True]};
sol = NDSolve[eq, f, {x, -10, 10}, {y, -10, 10}, {t, 0, 1}]

fu[x_, y_] = Abs@f[x, y, 1] /. sol;
Plot3D[fu[x, y], {x, -10, 10}, {y, -10, 10}, PlotRange -> All]