% Two-element truss, truss4,  3 nodes
clear                       % variables
format short e
% format long               % double precision

dtor = pi/180;              % convert degrees to radians
force = 50e3;               % applied force
fang = -60*dtor;            % angle of force
px = (force*cos(fang))
py = (force*sin(fang))
ang1 = 45
ang2 = 180-30
c1 = cos(ang1*dtor)
s1 = sin(ang1*dtor)
c2 = cos(ang2*dtor)
s2 = sin(ang2*dtor)

ex = [ 11e6 15e6 ]          % elastic modulus
area = [ 3 5 ]
len = [ 10 10*sqrt(2) ]
k12 = area(1)*ex(1)/len(1)  % element 1
k12c = k12*c1^2
k12s = k12*s1^2
k12cs = k12*c1*s1
k23 = area(2)*ex(2)/len(2)  % element 2
k23c = k23*c2^2
k23s = k23*s2^2
k23cs = k23*c2*s2
k33 = k12c + k23c
k44 = k12s + k23s
k34 = k12cs + k23cs
beta = k12 * 1e7
stiff = [ k12c   k12cs  -k12c   -k12cs    0      0  % stiff matrix
          k12cs  k12s   -k12cs  -k12s     0      0
         -k12c  -k12cs   k33     k34    -k23c  -k23cs
         -k12cs -k12s    k34     k44    -k23cs -k23s
          0        0    -k23c   -k23cs   k23c   k23cs
          0        0    -k23cs  -k23s    k23cs  k23s  ]
df =   [ 0; 0; px; py; 0; 0 ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(2,2)=zienk(2,2)+beta;
zienk(5,5)=zienk(5,5)+beta;
zienk(6,6)=zienk(6,6)+beta
format short e
displ = inv(zienk) * df;     % displacements
force = stiff * displ;       % force vector

% print the results
fprintf('\n    Two-element truss 6x6, truss4')
fprintf('\n    Your name, MENG 421, ')
disp(date)
fprintf('\n  Node    X Displ       FX       Y Displ       FY ')
fprintf('\n           inches       Kip       inches       Kip \n')
for i = 1:2:length(df)
  disp(sprintf('    %2.0f   %10.2e %7.1f    %11.2e %7.1f',...
   (i+1)/2,displ(i),force(i)/1000,displ(i+1),force(i+1)/1000))
end

% End of Matlab code

%    Two-element truss 6x6, truss4
%    Your name, MENG 421, 27-Feb-2004

%  Node    X Displ       FX       Y Displ       FY
%           inches       Kip       inches       Kip
%     1   -5.55e-010    18.3     -5.55e-010    18.3
%     2    2.84e-003    25.0     -1.39e-002   -43.3
%     3    1.31e-009   -43.3     -7.58e-010    25.0