Solving a time-independent 2D Schrödinger equation using PDE Toolbox

I want to solve the following PDE eigenvalue problem,

inside a square with Dirichlet boundary conditions on its edges, say for the first 250 eigenvalues and eigenfunctions. I don't want to implement this problem using the finite difference method, but rather prefer to use the already available PDE toolbox in MATLAB. When I looked at the examples presented in the MathWorks documentations for the PDE toolbox, I couldn't understand how to implement this problem using the PDE toolbox.

My question is; how can I implement this problem in MATLAB using the PDE toolbox? After that how can I access the eigenfunctions found by the code to take certain integrals involving these eigenfunction?