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


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.

  • 1
    $\begingroup$ EN doesn't seem to be defined - I assume that's just a typo here if it's working for you. $\endgroup$
    – N.J.Evans
    Commented Mar 23, 2023 at 13:19
  • $\begingroup$ Thanks, yes, corrected the typo $\endgroup$
    – Rom1984
    Commented Mar 23, 2023 at 13:28

2 Answers 2


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

 (*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]], 
 {kx, -1, 3}, {ky, -1, 3},
 {min, {-1, -1}, {4, 4}},
 {max, {2, 2}, {5, 5}}]

enter image description here enter image description here

  • $\begingroup$ 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 $\endgroup$
    – Rom1984
    Commented 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}]

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.