%% PROBLEM 3
clear
% set number of elements
N=5;
%define elemental properties  (A=pi*D^2/4)
EA=pi/4*[4.00^2;5.20^2;5.00^2;(5.00^2+4.00^2)/2;(4.00^2+3.00^2)/2];
L =[5.00;1.00;1.00;5.00;5.00];
%set up forcing vector
F =[1;0;0;0;0;-3];
%define boundary conditions
BC=[0;1;0;0;0;0];

%% SOLVER
%set up stiffness matrix
k=EA./L;
K=zeros(N+1,N+1);
for ii=1:N
K(ii,ii)=K(ii,ii)+k(ii);
K(ii+1,ii)=K(ii+1,ii)-k(ii);
K(ii,ii+1)=K(ii,ii+1)-k(ii);
K(ii+1,ii+1)=K(ii+1,ii+1)+k(ii);
end
%enforce BC's on [K] and {F}
K_reduced = K;
F_reduced = F;
for jj=N+1:-1:1
if BC(jj)==1
K_reduced(jj,:)=[];
K_reduced(:,jj)=[];
F_reduced(jj)=[];
end
end
%solve reduced system of equations
u_reduced=inv(K_reduced)*F_reduced;
%determine reaction at supports and plot results
r=1;
for jj=1:N+1
   if BC(jj)==1
       nodal_u(jj,1) = 0;
   else
       nodal_u(jj,1) = u_reduced(r);
       r=r+1;
   end
end
disp('Nodal displacement '); disp(nodal_u');
disp('Nodal force'); disp((K*nodal_u)');
disp('Stiffness matrix'); disp(K);