Why am I receiving this error?

span = 2000; 
numElements = 10; 
lengthPerElement = span / numElements;
pointLoad = 30 * 10^3; 
EI = 8 * 10^12; 

numNodes = numElements + 1; % Number of nodes (including the fixed end)
K_global = zeros(numNodes, numNodes)
F_load = zeros(numNodes, 1);

for i = 1:numElements
    ke = elementStiffness(EI, lengthPerElement);
    
    nodes = [i, i+1];
    K_global(nodes, nodes) = K_global(nodes, nodes) + ke;
    
    if i == numElements
        F_load(numNodes) = pointLoad;
    end
end

K_global(1, :) = 0;
K_global(1, 1) = 1;
F_load(1) = 0;

displacements = K_global \ F_load;

slope_at_free_end = (displacements(numNodes) - displacements(numNodes-1)) / lengthPerElement;

x_coordinates = linspace(0, span, numNodes);

figure;
plot(x_coordinates, displacements);
xlabel('X-coordinate (mm)');
ylabel('Deflection (mm)');
title('Deflection of Cantilever Beam');
grid on;

fprintf('Slope at the free end: %.4f radians\n', slope_at_free_end);


function ke = elementStiffness(EI, L)
ke = (EI / L^3) * [12, 6 * L, -12, 6 * L;
    6 * L, 4 * L^2, -6 * L, 2 * L^2;
    -12, -6 * L, 12, -6 * L;
    6 * L, 2 * L^2, -6 * L, 4 * L^2]
end

[Dyuman Joshi]

(My previous comment got deleted by mistake)

You are trying to add arrays of incompatible dimensions -

nodes = [i, i+1];
K_global(nodes, nodes) = K_global(nodes, nodes) + ke;

K_global(nodes,nodes) is a 2x2 and ke is 4x4.