% Set all constants
nx = 4;
ny = 4;
Lx = 1;
Ly = 1;
hbar = 1;
m = 1;
kx = @(nx) nx*pi/Lx;
ky = @(ny) ny*pi/Ly;
w = hbar*(kx(nx)^2+ky(ny)^2)/2*m;
A = sqrt(4/Lx*Ly);
[x,y] = meshgrid(0:Lx/100:Lx, 0:Ly/100:Ly);
t = linspace(0,10,500);

% Define the position and time functions
psi = @(x,y,nx,ny) A.*sin(kx(nx)*x).*sin(ky(ny)*y);
phi = @(t) exp(-1i.*w.*t);

filename = '2DParticle44.gif';
for j = 1:4:length(t)
    surf(x,y,real(psi(x,y,nx,ny)*phi(t(j))))
    xlim([0 Lx]), ylim([0 Ly]), zlim([-A A]);
    title('$\psi(x,y)$ for $(n_x,n_y) = (4,4)$','interpreter','latex')
    drawnow
    frame = getframe(1);
    im = frame2im(frame);
    [Q,map] = rgb2ind(im,256);
    if j == 1;
        imwrite(Q,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
    else
        imwrite(Q,map,filename,'gif','WriteMode','append','DelayTime',0);
    end
end