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$:[![enter image description here](https://i.stack.imgur.com/HYzUy.png)](https://i.stack.imgur.com/HYzUy.png)

with the following colormap

[![enter image description here](https://i.stack.imgur.com/77OqD.png)](https://i.stack.imgur.com/77OqD.png)

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

[![enter image description here](https://i.stack.imgur.com/4teBG.png)](https://i.stack.imgur.com/4teBG.png)

Answer: 
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", "|Ψ|"}]

```

[![enter image description here](https://i.stack.imgur.com/Y7m3P.png)](https://i.stack.imgur.com/Y7m3P.png)

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>).

[![enter image description here](https://i.stack.imgur.com/Y7m3P.png)](https://i.stack.imgur.com/Y7m3P.png)

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}}

```
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]

```

[![enter image description here](https://i.stack.imgur.com/wQLia.png)](https://i.stack.imgur.com/wQLia.png)