help correct the error

kindly help,

the error message reads 'index exceeds matrix dimensions

% Parameters
Lx = 1; % Length of domain in x-direction
Ly = 1; % Length of domain in y-direction
Lz = 1; % Length of domain in z-direction
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
T = 1; % Final time
alpha = 0.1; % Thermal diffusivity
dt = 0.001; % Time step size

% Grid spacings
dx = Lx/Nx;
dy = Ly/Ny;
dz = Lz/Nz;

% Create grid
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);

% Initial condition
u0 = sin(pi*X).*sin(pi*Y).*sin(pi*Z);

% Initialize solution
u = u0;

% Time-stepping loop
t = 0;
while t < T
    % Compute Laplacian using finite differences
    
    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
    d2u_dy2 = (-u(3:Nx, 2:Ny+2, 2:Nz) + 16*u(2:Nx, 2:Ny+1, 2:Nz) - 30*u(1:Nx, 2:Ny, 2:Nz) + 16*u(4:Nx, 2:Ny-1, 2:Nz) - u(5:Nx, 2:Ny-2, 2:Nz)) / 12*dy^2;
    d2u_dz2 = (-u(3:Nx, 2:Ny, 2:Nz+2) + 16*u(2:Nx, 2:Ny, 2:Nz+1) - 30*u(1:Nx, 2:Ny, 2:Nz) + 16*u(4:Nx, 2:Ny, 2:Nz-1) - u(5:Nx, 2:Ny, 2:Nz-2)) / 12*dz^2;
    
    % Update solution using the method of lines
    u(2:Nx, 2:Ny, 2:Nz) = u(2:Nx, 2:Ny, 2:Nz) + alpha*dt*(d2u_dx2 + d2u_dy2 + d2u_dz2);
    
    % Update time
    t = t + dt;
end

Index in position 1 exceeds array bounds. Index must not exceed 51.

% Plot final solution
figure;
slice(X, Y, Z, u, Lx/2, Ly/2, Lz/2);
xlabel('x');
ylabel('y');
zlabel('z');
colorbar;

Answer 1

[Walter Roberson]

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.