FEM Beam problem

% This Matlab code help you to calculate the displacements and reactions of
% Beam. You may Need to change the boundary conditions for different cases
% of Beam. If you are changing the number of elements, you may need to
% change the force vectors (F_udl & F_pl) in line 16 and 17.
clc; clear; close all;

L = 1; % Length in m
E = 2.1*10^8; % Modulus of Elasticity KN/m2
I = 2120/100^4; % Moment of Inertia m4
w = 1000; % Uniformly Distributed Load kN/m
wp = 1000; % Concentrated Load kN

BC = [1 2 5]; % Boundary Condition

F_udl = w*[-L/2;(L^2)/48-(L^2/48);(L^2/48)];
F_pl = wp*[-1;-0.05;0];

Ne = input('Type Number of element or press "Enter" to use default 2 Element ');
if (isempty(Ne))
Ne = 2;
else
end

Ndof = 2*(Ne+1);
Dof = 1:1:Ndof;

le = L/Ne;

% Compute the stiffness matrix for the structure
Kj = zeros(4,4);
Kj(1,1) = 6;
Kj(1,2) = 3*le;
Kj(1,3) = -6;
Kj(1,4) = 3*le;
Kj(2,2) = 2*le^2;
Kj(2,3) = -3*le;
Kj(2,4) = le^2;
Kj(3,3) = 6;
Kj(3,4) = -3*le;
Kj(4,4) = 2*le^2;

for ii = 2:4,
for jj = 1:ii,
Kj(ii,jj) = Kj(jj,ii);
end
end
Ke = (2*E*I/le^3).*Kj;

KG = zeros(Ndof,Ndof);
for ii = 1:Ne
aa = 2*ii-1;
KG(aa:aa+3,aa:aa+3) = KG(aa:aa+3,aa:aa+3)+Ke;
end

K_temp = KG;
for jj = 1:length(BC)
K_temp(:,BC(jj)) = 0;
K_temp(BC(jj),:) = 0;
end
K_temp(~any(K_temp,2),:) = [];
K_temp(:,~any(K_temp,1)) = [];
KR = K_temp;

FR = F_udl+F_pl;

UR = KR\FR;
UG = [0 0 UR(1) UR(2) 0 UR(3)]';
Rc = KG*(UG);

clc
disp('Displacement in mm & rad')
UR = [UR(1)*1000;UR(2);UR(3)];
disp(UR)
disp('Reaction in kN')
Reaction = [Rc(1)+w*L/4; Rc(2)+w*(L/2)^2/12; Rc(5)+w*L/4];
disp(Reaction)