function Beam_Ex
clc; clear all
% Beam example
F = 18; q = 10; EI = (210*10^6)*4*10^(-4);
L=2;
nodes =[0:L:3*L];
n=2*length(nodes);
debc=[1,2,7]; %Boundary conditions locations
ebcVals=zeros(length(debc),1); %Boundary conditions values
K=zeros(n); 
R = zeros(n,1); %Load Matrix
R(3)=-F;
% Generate equations for each element and assemble them.
for i=1:2
    lm=[2*(i-1)+1,2*(i-1)+2,2*(i-1)+3,2*(i-1)+4];
    [ke, rq] = BeamElement(2*EI, 0, nodes([i:i+1]));
    K(lm, lm) = K(lm, lm) + ke;
    R(lm) = R(lm) + rq;
end
for i=3
    lm=[2*(i-1)+1,2*(i-1)+2,2*(i-1)+3,2*(i-1)+4];
    [ke, rq] = BeamElement(EI, -q, nodes([i:i+1]));
    K(lm, lm) = K(lm, lm) + ke;
    R(lm) = R(lm) + rq;
end
K
R
% Nodal solution and reactions
d = NodalSoln(K, R, debc, ebcVals) %Displacement
va=[];
bma=[];
Va=[];
for i=1:2
    lm=[2*(i-1)+1,2*(i-1)+2,2*(i-1)+3,2*(i-1)+4];
    [v, bm, V]=BeamResults(2*EI, 0, nodes([i:i+1]), d(lm));
    va = [va; v]; bma = [bma; bm]; Va = [Va; V];
end
for i=3
    lm=[2*(i-1)+1,2*(i-1)+2,2*(i-1)+3,2*(i-1)+4];
    [v, bm, V]=BeamResults(EI, -q, nodes([i:i+1]), d(lm));
    va = [va; v]; bma = [bma; bm]; Va = [Va; V];
end
va
bma
Va
node_plot= va(:,1); 

 
node_plot
va(:,2)
figure (1)
plot(node_plot,va(:,2),'o-','LineWidth',2,'MarkerSize',10),title('Verical displacement'), xlabel('x'),ylabel('v')
grid on
figure (2)
plot(node_plot,bma(:,2),'o-','LineWidth',2,'MarkerSize',10),title('Moment Diagram'), xlabel('x'),ylabel('M')
grid on
figure (3)
plot(node_plot,Va(:,2),'o-','LineWidth',2,'MarkerSize',10),title('Shear Diagram'), xlabel('x'),ylabel('M')
grid on

function [ke, rq] = BeamElement(EI, q, coord)
% [ke, rq] = BeamElement(EI, q, coord)
% Generates equations for a beam element
% EI = beam stiffness
% q = distributed load
% coord = coordinates at the element ends

L=coord(2)-coord(1);
ke = [(12*EI)/L^3, (6*EI)/L^2, -((12*EI)/L^3), (6*EI)/L^2;
    (6*EI)/L^2, (4*EI)/L, -((6*EI)/L^2), (2*EI)/L;
    -((12*EI)/L^3), -((6*EI)/L^2), (12*EI)/L^3, -((6*EI)/L^2);
    (6*EI)/L^2, (2*EI)/L, -((6*EI)/L^2), (4*EI)/L];
rq = [(L*q)/2; (L^2*q)/12; (L*q)/2; -((L^2*q)/12)];

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

function [v, bm, V] = BeamResults(EI, q, coord, dn)
% [v, bm, V] = BeamResults(EI, q, coord, dn)
% Generates beam element results
% EI = beam stiffness
% q = distributed load
% coord = coordinates at the element ends
% dn = nodal solution
L=coord(2)-coord(1);
v=[]; bm=[]; V=[];
% Change increment to get results at more points
for s=0:0.01:L
    n = [(2*s^3)/L^3 - (3*s^2)/L^2 + 1, s^3/L^2 - (2*s^2)/L + s, ...
            (3*s^2)/L^2 - (2*s^3)/L^3, s^3/L^2 - s^2/L];
    v = [v; [coord(1)+s, n*dn+(q*(L - s)^2*s^2)/(24*EI)]];
    dn2 = [(12*s)/L^3 - 6/L^2, (6*s)/L^2 - 4/L, 6/L^2 - (12*s)/L^3, ...
            (6*s)/L^2 - 2/L];
    bm = [bm; [coord(1)+s, EI*dn2*dn+(q*(L^2 - 6*s*L + 6*s^2))/(12)]];
    dn3 = [12/L^3, 6/L^2, -(12/L^3), 6/L^2];
    V = [V; [coord(1)+s, EI*dn3*dn+((q*(12*s - 6*L))/(12))]];
end
