% Axial1, stepped bar, 2 elements, 3 nodes
%
clear                          % variables
format short e                 % E format
% format long                  % double precision
p1 = 1000                      % load at node 1
% define elastic modulus, length, and diameter
% the following 3 vectors need one entry per element
ex = [ 29 11 ]*1e6             % elastic modulus
len = [10 12]                  % element lengths
dia = [1.1  1.6]               % element diameters
a = pi/4*dia .^2               % area

% define the stiffness matrix for each element
k12 = a(1)*ex(1)/len(1)        % stiff const
k23 = a(2)*ex(2)/len(2)
k22 = k12 + k23
beta = k12*1e7                 % Zienkiewicz factor

stiff = [ k12    -k12     0    % stiffness matrix
         -k12     k22   -k23
          0      -k23    k23 ]
df =    [ p1;  0; 0*beta ]     % combined d, F
zienk = stiff;
zienk(3,3)=zienk(3,3)+beta     % Zienkiewicz method

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

for i = 1:length(df)-1         % element displacements
  eldis(i) = displ(i+1) - displ(i);
end
stress = ex .* eldis ./ len;   % element stresses

% print the results
fprintf('\n           Axial1 ')
fprintf('\n  Your name, MENG 421, ')
disp(date)

fprintf('\n   Node Displacement    Force')
fprintf('\n            Inches       Lb \n')
for i = 1:length(df)
  fprintf('    %2.0f   %11.3e %7.0f \n',i,displ(i),force(i))
end

fprintf('\n   Elem    Stress')
fprintf('\n             Psi \n')
for i = 1:length(df)-1
  fprintf('    %2.0f   %7.0f \n',i,stress(i))
end

% End of Matlab code

% The following is the output from this program

%          Axial1
%   Your name, MENG 421, 05-Feb-2004

%   Node Displacement    Force
%            Inches       Lb
%     1    9.054e-004    1000
%     2    5.426e-004       0
%     3    3.628e-011   -1000

%   Elem    Stress
%             Psi
%     1     -1052
%     2      -497