You could use the principal stress function template below to modify the call to
eig(...) and return the eigenvectors too.
function [ se_p ] = princse( se_x, se_y, se_z, se_xy, se_yz, se_xz, icomp )
% PRINCSE Compute principal stresses and strains.
%
% [ SE_P ] = PRINCSE( SE_X, SE_Y, SE_Z, SE_XY, SE_YZ, SE_XZ, ICOMP ) Computes
% principal stresses and strains (for the ICOMP component, default 1). In two
% dimensions only inputs SE_X, SE_Y, SE_Z, SE_XY, and ICOMP are required.
n = length( se_x );
eg = zeros( 3, 3, n );
eg(1,1,:) = se_x;
eg(2,2,:) = se_y;
eg(3,3,:) = se_z;
eg(2,1,:) = se_xy;
eg(1,2,:) = se_xy;
if( nargin>=6 )
eg(3,2,:) = se_yz;
eg(2,3,:) = se_yz;
eg(3,1,:) = se_xz;
eg(1,3,:) = se_xz;
if( nargin<7 )
icomp = 1;
end
else
if( nargin<5 )
icomp = 1;
else
icomp = se_yz;
end
end
d = zeros( 3, n );
for i=1:n
d(:,i) = eig( eg(:,:,i) );
end
srt = sort( [ d(1,:); d(2,:); d(3,:) ] );
srt = srt( [3;2;1], : );
se_p = reshape( srt(icomp,:), size(se_x) );