Nx = 50; % Number of grid points in x-direction
Ny = 50; % Number of grid points in y-direction
Nz = 50; % Number of grid points in z-direction
50's.
x = linspace(0, Lx, Nx+1);
y = linspace(0, Ly, Ny+1);
z = linspace(0, Lz, Nz+1);
[X, Y, Z] = meshgrid(x, y, z);
50+1 by 50+1 by 50+1 --> 51 x 51 x 51
d2u_dx2 = (-u(3:Nx+2, 2:Ny, 2:Nz) + 16*u(2:Nx+1, 2:Ny, 2:Nz) - 30*u(1:Nx, 2:Ny, 2:Nz) + 16*u(4:Nx-1, 2:Ny, 2:Nz) - u(5:Nx-2, 2:Ny, 2:Nz)) /12*dx^2; %HERE
Nx+2 --> 50+2 --> 52. But the array is only 51.