%% PROBLEM 2c
clc
clear
for num_elements=3:500
% set number of elements
N=num_elements;
%define elemental properties
EA=[1  ;1  ;(2+([0:N-3]'+[1:N-2]')/(N-2))/2];
L =[0.5;0.5;ones(N-2,1)/(N-2)];
%set up forcing vector
F =[1;0;zeros(N-1,1)];
%define boundary conditions
BC=[0;0;zeros(N-2,1);1];

%% 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
displ_first_node(num_elements,1) = num_elements;
displ_first_node(num_elements,2) = nodal_u(1);
end
disp('Num Elements | Displ. Node 1'); disp(displ_first_node(3:num_elements,:));
