% Simply supported beam, two nodes, simpm.m
% W8x31 beam
delete simpm.d             % remove old file
clear                      % variables
format short e
m2 = -250e3                % moment at free end 

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

stiff = [ L12    L6   -L12     L6     % stiff matrix
          L6     L4   -L6      L2
         -L12   -L6    L12    -L6
          L6     L2   -L6      L4 ]
df =   [ 0*beta; 0; 0*beta; m2 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(3,3)=zienk(3,3)+beta

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

% print the results
fprintf('\n   Simply supported beam, simpm')
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

%   Simply supported beam, simpm
%   Your name, MENG 421, 20-Apr-2004

%   Node   Deflection    Rotation      Force    Moment
%     1     6.31e-009     1.26e-003   -2500     5.82e-011 
%     2    -6.31e-009    -2.53e-003    2500    -2.50e+005