% Simply supported beam, three nodes, simp2.m
% W8x31 beam, conc force, 6 kips
clear                      % variables
format short e
f2 = -6e3                  % force at center

ex = 30e6                  % steel elastic modulus
ei = ex*110; s = 27.5      % E x I and section modulus

Len = 100
L2 = 2*ei/Len; L4 = 4*ei/Len
L6 = 6*ei/Len^2; L12 = 12*ei/Len^3
beta = L12 * 1e7

stiff1 = [ L12    L6   -L12     L6     % stiff matrix
          L6     L4   -L6      L2
         -L12   -L6    L12    -L6
          L6     L2   -L6      L4 ]

stiff = zeros(6); stiffa = zeros(6);
stiff(1:4,1:4)= stiff1;
stiffa(3:6,3:6)= stiff1;
stiff = stiff + stiffa
df =   [ 0*beta; 0; f2; 0; 0*beta; 0 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(5,5)=zienk(5,5)+beta

format short e
displ = inv(zienk) * df;     % displacements
displt = displ'
force = (stiff * displ)'     % force vector

% print the results
fprintf('\n   Simp2, simply-supported beam, three nodes  ')
fprintf('\n      Your name, MENG 421, ')
disp(date)
fprintf('\n   Node    Deflec     Rotation      Force    Moment\n')
for i = 1:2:length(df-2)
  disp(sprintf('    %2.0f   %8.4f   %11.2e   %6.0f   %8.0f',...
    (i+1)/2,displ(i),displ(i+1),force(i),force(i+1)))
end

%   Simp2, simply-supported beam, three nodes
%      Your name, MENG 421, 05-Feb-2004

%   Node    Deflec     Rotation      Force    Moment
%     1    -0.0000    -4.55e-003     3000          0
%     2    -0.3030    -2.85e-019    -6000          0
%     3    -0.0000     4.55e-003     3000          0