finite element to calculate shear force diagram and bending moment diagram deflection of cantilever and distributed load in matlab

main.m
% Define properties
E = 2.1e11; % Young's modulus (Pa)
I = 4.0833e-6; % Moment of inertia (m^4)

% Element properties
L = 3; % Length of the beam (m)
numElements = 10; % Number of elements
elementLength = L / numElements;

% Distributed load
w = 1000; % Distributed load (N/m)

% Initialize global matrices
K = zeros(numElements + 1, numElements + 1);
F = zeros(numElements + 1, 1);

% Assemble stiffness matrix and load vector
for i = 1:numElements
    ke = (E*I/elementLength^3) * [12, 6*elementLength, -12, 6*elementLength; 
                                  6*elementLength, 4*elementLength^2, -6*elementLength, 2*elementLength^2; 
                                  -12, -6*elementLength, 12, -6*elementLength; 
                                  6*elementLength, 2*elementLength^2, -6*elementLength, 4*elementLength^2];
    fe = [w*elementLength/2; w*elementLength^2/12; w*elementLength/2; -w*elementLength^2/12];
    
    % Assemble into global matrices
    K(i:i+1, i:i+1) = K(i:i+1, i:i+1) + ke;
    F(i:i+1) = F(i:i+1) + fe;
end

% Apply boundary conditions
K(1,:) = 0;
K(1,1) = 1;
F(1) = 0;

% Solve for displacements
U = K \ F;

% Calculate shear force diagram, bending moment diagram, and deflection
V = diff(U) .* (E*I/elementLength^3);
M = -E*I*diff(diff(U)) / elementLength^2;
D = U(2:end);

% Plot the diagrams
x = linspace(0, L, numElements+1);
figure;
subplot(3,1,1); plot(x, V); title('Shear Force Diagram');
subplot(3,1,2); plot(x, M); title('Bending Moment Diagram');
subplot(3,1,3); plot(x, D); title('Deflection');

1538 chars
50 lines