% Simply supported beam, three nodes, simpd2.m
% W8x31 beam, uniform load = 60 lb/in
clear                       % variables
w2 = 60                     % uniform load

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

Len = 100;
f2 = w2*Len; m3 = f2 * Len / 12;
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; -m3; -f2; 0; 0*beta; m3 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(5,5)=zienk(5,5)+beta

displ = inv(zienk) * df;       % displacements
force = stiff * displ;         % force vector

% print the results
fprintf('\n  Simpd2, simply-supported beam, 3 nodes ')
fprintf('\n      Your name, MENG 421, ')
disp(date)
fprintf('\n   Node  Deflection  Rotation     Force    Moment')
fprintf('\n          Inches      Radians     Kips     Kip-in \n')
for i = 1:2:length(df-2)
  fprintf('    %2.0f   %6.3f   %11.2e   %7.3f   %6.2f \n',...
    (i+1)/2,displ(i),displ(i+1),force(i)/1000,force(i+1)/1000)
end

%   Simpd2, simply-supported beam, 3 nodes
%      Your name, MENG 421, 20-Feb-2004

%   Node  Deflection  Rotation     Force    Moment
%          Inches      Radians     Kips     Kip-in
%     1   -0.000    -6.06e-003     3.000    -50.0
%     2   -0.379    -3.52e-019    -6.000     -0.0
%     3   -0.000     6.06e-003     3.000     50.0