% Axial1, stepped bar, 2 elements, 3 nodes, both ends fixed
clear                          % variables
format short e                 % E format
% format long                  % double precision

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 =    [ 0*beta; 1000; 0*beta] % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;    % Zienkiewicz method
zienk(3,3)=zienk(3,3)+beta

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           Axial2 ')
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

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

%   Node Displacement    Force
%            inches       Lb
%     1    2.174e-011    -599
%     2    2.174e-004    1000
%     3    1.454e-011    -401

%   Elem    Stress
%             Psi
%     1       631
%     2      -199