How to solve analytically a 2D Sturm-Liouville Problem?

Question

Maybe it is a syntax problem, since I haven't been using Mathematica for much time, but I haven't been able to find a similar example in internet.

I am just trying to solve analitically the Sturm-Liouville problem of Schrödinger equation for a particle in a 2D box $[0,1 ] \times [0,1]$. Id est:

$ - \frac{\hbar^2}{2m} \Delta \psi(x,y) = E \psi(x, y) $

$\psi(0,y) = \psi(1,y) = \psi(x,0) = \psi(x,1) = 0$

This is what I wrote:

eqn = 2* Pi^2 *z[x, y] + (D[z[x, y], {x, 2}] + D[z[x, y], {y, 2}]) == 0;
bc = {z[x, 0] == 0, z[x, 1] == 0, z[0, y] == 0, z[1, y] == 0};
sol = DSolve[{eqn, bc}, z[x, y], {x, y}]

And this is what I get:

Answer 1

You can do this with M, but you need to help it. Start with the pde

 pde = D[ψ[x, y], x, x] + D[ψ[x, y], y, y] == -2 m En ψ[x, y]/ℏ^2

Separate variables

ψ[x_, y_] = X[x] Y[y]

pde = pde/ψ[x, y] // Expand

D[X[x], x, x]/X[x] + D[Y[y], y, y]/Y[y] == -((2 En m)/ℏ^2)

The first term is a fn of x only, the second a function of y and the rhs is a constant, so each term must be equal to a constant. A negative constant will give the desired sinusoidal solutions.

DSolve[{pde[[1, 1]] == -a^2, X[0] == 0}, X[x], x] // Flatten

{X[x] -> C[2] Sin[a x]}

Only put in one condition because, if not, it will return only the trivial solution of 0. To solve for the condition at x = 1 for a

X[x_] = X[x] /. % /. C[2] -> c1

Reduce[X[1] == 0, a]

    (C[1] ∈ Integers && (a == 2 π C[1] || a == 2 π C[1] + π)) || c1 == 0

which boils down to

a = n Pi

and

$Assumptions = n ∈ Integers

Now the y part. Subst the x part into the pde.

pde

D[Y[x], y, y]/Y[y] - n^2 π^2 == -((2 En m)/ℏ^2)

Set the Y term to another negative constant.

DSolve[{D[Y[y], y, y]/Y[y] == -b^2, Y[0] == 0}, Y[y], y] // Flatten

{Y[y] -> C[2] Sin[b y]}

Y[y_] = Y[y] /. % /. C[2] -> 1

Sin[b y]

I set the C[2] to 1 because it is multiplied by the constant in the x equation leaving only one constant anyway. The condition at y = 1

Y[1] == 0

Sin[b] == 0

and as before

b = k π

$Assumptions = $Assumptions && k ∈ Integers

pde

-(π^2 k^2) - π^2 n^2 == -((2 En m)/ℏ^2)

Solve for the Energies

Solve[pde, En] // Flatten

{En -> (ℏ^2 (π^2 k^2 + π^2 n^2))/(2 m)}

Normalize the wave function

normint = Integrate[ψ[x, y]^2, {x, 0, 1}, {y, 0, 1}] == 1

c1^2/4 == 1

c1 = c1 /. Solve[normint, c1][[2]]

ψ[x, y]

2 Sin[π k y] Sin[π n x]

Answer 2

As mentioned in the comments, the equation and boundary conditions need to be in the same list.

Beyond that, Mathematica has a little bit of difficulty solving the exact equation you've given: but if you're happy to sacrifice a bit of exactitude in the name of getting the right answer, you can try

eqn = (1 + $MachineEpsilon)*2*Pi^2*z[x, y] + (D[z[x, y], {x, 2}] + D[z[x, y], {y, 2}]) == 0; 
bc = {z[x, 0] == 0, z[x, 1] == 0, z[0, y] == 0, z[1, y] == 0}; 
sol = z[x, y] /. DSolve[Join[{eqn}, bc], z[x, y], {x, 0, 1}, {y, 0, 1}]

which gives

{2*Inactive[Sum][0, {K[1], 1, Infinity}]}

i.e., 0.