a code using function_handle is not woring

matlab

% Define the domain and grid size 
L = 3; % length of the domain 
N = 2; % number of grid points in each direction 
x = linspace(-L/2,L/2,N); 
y = linspace(-L/2,L/2,N); 
z = linspace(-L/2,L/2,N); 
[X,Y,Z] = meshgrid(x,y,z); 
dx = x(2)-x(1); 

% Define the potential function 
V = -1./sqrt(X.^2+Y.^2+Z.^2); % Coulomb potential 

% Define the Laplacian operator using central differences 
Lap = @(f) (circshift(f,[0,0,-1])+circshift(f,[0,0,1])+circshift(f,[0,-1,0])+circshift(f,[0,1,0])+circshift(f,[-1,0,0])+circshift(f,[1,0,0])-6*f)/(dx^2)*(-0.5); 
% Set up the initial wave function (1s state) 
psi = exp(-sqrt(X.^2+Y.^2+Z.^2)); % Gaussian wave packet 

% Normalize the wave function 
psi = psi./sqrt(sum(abs(psi(:)).^2*dx^3)); 

% Define the Hamiltonian operator 
H=Lap+diag(V(:));

Operator '+' is not supported for operands of type 'function_handle'.

% Set up the time evolution parameters 
dt = 0.01; % time step 
%tmax = 10; % maximum time 
tmax = 0.02; % maximum time
t = 0:dt:tmax; 

% Solve the time-dependent Schrodinger equation using the Crank-Nicolson method 
for n = 1:length(t) 
    psi = (eye(N^3)-0.5i*dt*H)\(eye(N^3)+0.5i*dt*H)*psi; 
    psi = psi./sqrt(sum(abs(psi(:)).^2*dx^3)); % renormalize 
end 

% Calculate the energy levels by finding the eigenvalues of the Hamiltonian 
[E,D] = eigs(H,10,'sr'); 
E = diag(E); % extract the eigenvalues 
E = E(E<0); % take only the negative energy levels (bound states) 

% Print the energy levels 
disp('Energy levels:'); 
disp(E); 

[Torsten]

"Lap" is a function handle and has to be called with an explicit argument f:

H=Lap(f)+diag(V(:));

I don't know what f is in your code from above.

[Stephen23]

Duplicate:

https://www.mathworks.com/matlabcentral/answers/2165562-operator-is-not-supported-for-operands-of-type-function_handle