% Simply supported beam, two nodes, simpd.m
% W8x31 beam, uniform load = 120 lb/in
clear                      % variables
format short e
w2 = 120                   % uniform load

ex = 30e6                  % steel elastic modulus
ei = ex*110                % E x I
l = 100;    f2 = w2*l/2; m2 = f2 * l / 6;
l2 = 2*ei/l; l4 = 4*ei/l
l6 = 6*ei/l^2; l12 = 12*ei/l^3
beta = l12 * 1e7
stiff = [ l12    l6   -l12     l6     % stiff matrix
          l6     l4   -l6      l2
         -l12   -l6    l12    -l6
          l6     l2   -l6      l4 ]
df =   [ 0*beta; -m2; 0*beta; m2 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(3,3)=zienk(3,3)+beta
format short e
displ = inv(zienk) * df;     % displacements
displt = displ'
force = (stiff * displ)'     % force vector

fprintf('\n  Simpd, simply-supported beam, two nodes, 4x4 ')
fprintf('\n      Your name, MENG 421, ')
disp(date)
fprintf('\n   Node   Deflection    Rotation      Force    Moment\n')
for i = 1:2:length(df-2)
  fprintf('    %2.0f   %11.2e   %11.2e   %5.0f   %11.2e \n',...
    (i+1)/2,displ(i),displ(i+1),force(i),force(i+1))
end

%    Simpd, simply-supported beam, two nodes, 4x4
%      Your name, MENG 421, 05-Feb-2004

%   Node   Deflection    Rotation      Force    Moment
%     1    -4.14e-025    -1.52e-003       0    -1.00e+005
%     2     0.00e+000     1.52e-003       0     1.00e+005