% Cantilever beam, two nodes, cant1.m
% W8x31 beam
clear                      % variables
format short e
f2 = -3000                 % force 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*beta; f2; 0 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(2,2)=zienk(2,2)+beta;

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

% print the results
fprintf('\n   Cantilever beam, two nodes, cant1')
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   %6.2f   %5.1f \n',...
    (i+1)/2,displ(i),displ(i+1),force(i)/1000,force(i+1)/1000)
end

%    Cantilever beam, two nodes, cant1
%    Your name, MENG 421, 16-Apr-2004

%    Node  Deflection  Rotation    Force  Moment
%           Inches     Radians      Kips   Kip-in 
%      1   -0.000    -7.58e-007     3.00   300.0 
%      2   -0.303    -4.55e-003    -3.00     0.0