Andrei's solution below works, although I had to make a couple of tweaks to get the behaviour just right so that it works in general for any 3x3xL multidimensional array:
%{
Construct a 3x3xL multidimensional array of the
stress tensor history, where L is the history
length
%}
% Direct stress components
normals = [Sxx; Syy; Szz]';
% Shear stress components
shears = [Txy; Txz; Tyz]';
% Assign normal stress component histories to diagonal terms
diagonals = repmat(eye([3.0, 3.0]), 1.0, 1.0, L);
diagonals(diagonals > 0.0) = normals';
% Assign shear stress component histories to non-diagonal terms
nonDiagonals = repmat(tril(ones(3.0), -1.0), 1.0, 1.0, L);
nonDiagonals(nonDiagonals > 0.0) = shears';
% Combine diagonal and non-diagonal terms
multidimensionalStressTensor = diagonals + nonDiagonals +
permute(nonDiagonals, [2.0, 1.0, 3.0]);
%{
Get the Eigenvalues for each frame of the tensor data. This code
was written by Bruno Luong and can be found on the exchange (File
ID: #27680)
%}
eigenvalues = eig3(multidimensionalStressTensor);
%{
The principal stresses are the maximum, median and
minimum values of the Eigenvector at each tensor
frame
%}
s1 = max(eigenvalues);
s2 = median(eigenvalues);
s3 = min(eigenvalues);
Initial tests indicate that this method is approximately 4 times faster than a FOR loop with 2D matrices.
