clc
clear 
close all 
%Data
FA=15;%kN
FC=80;%kN
FE=10;%kN
LAB=1;%m
LBC=1;%m
LCD=1;%m
LDE=1;%m
L=4;%m
E=200*10^6;%kPa
I=24*10^-6;%m^4
%SOLUTION
%Reaction
syms RB RD %Unknowns 
%Sum of moment at B
equ1=-FA*LAB+FC*LBC-RD*(LBC+LCD)+FE*(LBC+LCD+LDE)==0;
%Sum of moment at D
equ2=-FA*(LAB+LBC+LCD)+RB*(LBC+LCD)-FC*LCD+FE*LDE==0;
%Solving RA and RC
[RB RD]=solve([equ1 equ2],[RB RD]);
RB=double(RB)
RD=double(RD) %Range of x
dx=0.001; % Increment for discretization of beam length
x=0:dx:L; % range variable of beam length
%processing BM and SF programming loop at each section
for i=1:length(x)
 % section AB
if 0<=x(i)&&x(i)<LAB
 M(i)=-FA*x(i);
 V(i)=-FA;
 % Section BC
elseif LAB<=x(i)&&x(i)<LAB+LBC
 M(i)=-FA*x(i)+RB*(x(i)-LAB);
 V(i)=-FA+RB;
 % Section CD
elseif LAB+LBC<=x(i)&&x(i)<LAB+LBC+LCD
 M(i)=-FA*x(i)+RB*(x(i)-LAB)-FC*(LAB+LBC+LCD);
 V(i)=-FA+RB-FC;
 %Section DE
 M(i)=-FA*x(i)+RB*(x(i)-LAB)-FC*(LAB+LBC+LCD)+RD*(LAB+LBC+LCD+LDE);
 V(i)=-FA+RB-FC+RD;
end
end
figure(1)
subplot(2,1,1)
bar(x,V);
title('Shear Force Diagram')
xlabel('Length [m]')
ylabel('SF [kN]')
grid on
subplot(2,1,2)
bar(x,M);
title('Bending Moment Diagram')
xlabel('Length [m]')
ylabel('BM [kNm]')
grid on;
%Maximum shear force
Vmax=max(V)
%Maximum bending moment
Mmax=max(M)
%Numerical Integration
slope0=cumtrapz(M)*dx/(E*I);
deflect0=cumtrapz(slope0)*dx;
%Boundary conditions
%n=((last-first)/dx)+1
x1=0;
x2=4;
n1=((x1-0)/dx)+1;
n2=((x1-0)/dx)+1;
%Finding constant C
C=-[x1,1;x2,1]\[deflect0(n1);deflect0(n2)];
%Slope and deflection equation
slope=slope0+C(1);
y=deflect0+C(1)*x+C(2);
figure(2)
subplot(2,1,1)
plot(x,slope);
title('Slope Diagram')
xlabel('Length [m]')
ylabel('slope [rad]')
grid on
subplot(2,1,2)
plot(x,y);
title('Deflection Diagram')
xlabel('Length [m]')
ylabel('Deflection [m]')
grid on
xlim=ylim