how to construct the 2D finite-difference representation of Laplacian by Matlab

Question

I've been trying to solve the following Schrodinger equation numerically,

\begin{equation} \begin{split} \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} + U(x,y) \psi =E \psi \\ U(x,y) =cos(x)cos(y) \\ \end{split} \end{equation}

There is a Matlab codes (two-dimensional Schrodinger equation),

%***********************************************************************
% Program: Calculate two-dimensional time-independent schrodinger equation
%***********************************************************************
clear;
%% Parameters for solving problem in the interval 0 < x < L
L = pi/2;                   % Interval Length
N = 51;                    % No. of coordinate points
x = linspace(0,L,N)';       % Coordinate vector
dx = x(2) - x(1);           % Coordinate step

%% Construct the Potential matrix
    %v = @(x,y) 10*(sinh(x)^2 + sinh(y)^2)/(cosh(x) + cosh(y))^2;
     v = @(x,y) cos(x)*cos(y);
    U = zeros(N*N,N*N);
    for i = 1:N
       for j = 1:N
          pos = j + N*(i - 1);
          U(pos, pos) = v(-L/2 + dx*(i-1), -L/2 + dx*(j-1)); 
       end
    end



%% Five-point finite-difference representation of Laplacian
e   = ones(N,1);
spe = spdiags([e -2*e e], -1:1,N,N);
Iz  = speye(N);
Lap  = kron(Iz,spe)+kron(spe,Iz);
%full(Lap)
% 
% Next modify Lap so that it is consistent with f(0) = f(L) = 0
%Lap(1,1) = 0; Lap(1,2) = 0; Lap(2,1) = 0;    % So that f(0) = 0
%Lap(N,N-1) = 0; Lap(N-1,N) = 0; Lap(N,N) = 0;% So that f(L) = 0

% Total Hamiltonian
hbar = 1; m = 1;      % constants for Hamiltonian
H = -1/2*(hbar^2/m)*Lap + U;

%% Find lowest nmodes eigenvectors and eigenvalues of sparse matrix
[Psi, E] = eig(H);
    eval = diag(E);

    [useless, perm] = sort(eval);
    eval = eval(perm); Psi = Psi(:,perm);

  %% Plot a choice of Psi
    fprintf('Preparing the plot... \n');
    k = 0;
    PsiFn = zeros(N,N);

    for j = 1:N
        for i = 1:N
            PsiFn(i,j) = Psi(i + N*(j - 1), k + 1);
        end
    end
    [y,z] = meshgrid(-L/2:dx:L/2);
    % Ploting |Psi|^2
    figure;
   % surf(y,z,PsiFn);
surf(y,z,abs(PsiFn),'EdgeColor','none','LineStyle','none','FaceLighting','phong');

but this is only for N <51（The computer is out of memory. ） could anyone tell me a better way?