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

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:

• You must have the equation and the initial/boundary conditions in a list. This however does not mean that DSolve can solve this equation. That would be a separate issue ! – Lotus Oct 5 '17 at 12:22
• If you and @Miguel Bolín are the same person, please contact moderator for merging the accounts. – Alexey Popkov Oct 5 '17 at 12:43

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. You can do this with M, but you need to help it. Start with the pde  pde = D[\[Psi][x, y], x, x] + D[\[Psi][x, y], y, y] == -2 m En \[Psi][x, y]/\[HBar]^2  Separate variables \[Psi][x_, y_] = X[x] Y[y] pde = pde/\[Psi][x, y] // Expand D[X[x], x, x]/X[x] + D[Y[y], y, y]/Y[y] == -((2 En m)/\[HBar]^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] \[Element] Integers && (a == 2 \[Pi] C[1] || a == 2 \[Pi] C[1] + \[Pi])) || c1 == 0  which boils down to a = n Pi  and $Assumptions = n \[Element] Integers


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

pde

D[Y[x], y, y]/Y[y] - n^2 Pi^2 == -((2 En m)/\[HBar]^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 Pi

$Assumptions =$Assumptions && k \[Element] Integers

pde

-(\[Pi]^2 k^2) - \[Pi]^2 n^2 == -((2 En m)/\[HBar]^2)


Solve for the Energies

Solve[pde, En] // Flatten

{En -> (\[HBar]^2 (\[Pi]^2 k^2 + \[Pi]^2 n^2))/(2 m)}


Normalize the wave function

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

c1^2/4 == 1

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

\[Psi][x, y]

2 Sin[\[Pi] k y] Sin[\[Pi] n x]