# Setting up a PDE for a time-independent Schrödinger equation

I am trying to solve a 2d Schrödinger equation with a non-separable potential because I want to calculate the probability of reflection for different angles of incidence of a plane wave. Here I exclude the potential as I have a problem with setting up the basics (mostly boundary conditions). I have been reading the documentation and I am still a little bit confused about whether I need to use a NeumannValue.

My attempt so far, I managed to set boundary conditions for a wave that travels on the diagonal of the domain but only when the direction is non-normalised, otherwise this fails, and also doesn't work for any other direction.

(*Define the domain*)
xMin = 0;
xMax = 10;
yMin = -5;
yMax = 5;

(*Define the Schrödinger equation*)
EN = 2
eqn = {-Laplacian[u[x, y], {x, y}] == EN u[x, y]};

(*Define the plane wave boundary condition*)
k = {1, -1}

bc = {u[x, yMin] == Exp[I k . {x, yMin}],
u[x, yMax] == Exp[I k . {x, yMax}],
u[xMin, y] == Exp[-I k . {xMin, y}],
u[xMax, y] == Exp[-I k . {xMax, y}]};

(*Solve the Schrödinger equation with the boundary condition*)
sol = NDSolve[{eqn, bc}, u, {x, xMin, xMax}, {y, yMin, yMax}];

(*Plot the solution*)
Plot3D[Re[u[x, y] /. sol], {x, xMin, xMax}, {y, yMin, yMax}]


Question

How do I modify the boundary condition so this works for any arbitrary direction? Do I need to specify conditions for a reflection (presumably NeumannValue)? I assumed that the wave will reflect off of the potential so if the potential is not there once it reaches the wall it should go through.

• EN doesn't seem to be defined - I assume that's just a typo here if it's working for you. Mar 23, 2023 at 13:19
• Thanks, yes, corrected the typo Mar 23, 2023 at 13:28

Apparently, the code works fine for any of boundary conditions and domain:

nrg = 2;
eqn = {-Laplacian[u[x, y], {x, y}] == nrg u[x, y]};

Manipulate[
(*Define the plane wave boundary condition*)
k = {kx, ky};
bc = {
u[x, min[[2]]] == Exp[I k.{x, min[[2]]}],
u[x, max[[2]]] == Exp[I k.{x, max[[2]]}],
u[min[[1]], y] == Exp[-I k.{min[[1]], y}],
u[max[[1]], y] == Exp[-I k.{max[[1]], y}]};

(*Solve the Schrodinger equation with the boundary condition*)
sol = NDSolve[{eqn, bc},
u, {x, min[[1]], max[[1]]}, {y, min[[2]], max[[2]]}];

(*Plot the solution*)
Plot3D[Re[u[x, y] /. sol], {x, min[[1]], max[[1]]}, {y, min[[2]],
max[[2]]}],
{kx, -1, 3}, {ky, -1, 3},
{min, {-1, -1}, {4, 4}},
{max, {2, 2}, {5, 5}}]


• The solution behaves a bit weird once you make the domain larger. I have normalised the wave vector but once the x and y axis are increased to over 10, small artefacts begin to appear on the wave peaks, and when they're close to 30, the solution is no longer a plane wave. I wonder if this has something to do with the discretization or the method used Mar 25, 2023 at 15:59

The energy and k are not independent. From your code I take it that you set hbar=1 and 2m=1. Then the norm of k must be Sqrt[E].

Therefore for an arbitrary direction, we may choose k with arbitrary values of nx and ny:

k = Sqrt[2] Normalize[{nx,ny}]


E.g. for a wave along the x axis:

k = Sqrt[2] Normalize[{1, 0}]


• Great! Thanks very much! Yes, I have lost the forest for the trees, you're absolutely correct, it works fine now! Mar 23, 2023 at 14:33