function X = MatTriSol(S,T,B)
% S (m-by-m) and T (n-by-n) are upper triangular and B is an m-by-n.
% X is an m-by-n matrix such that SXT +  X = B.

% Compare first columns:  S*X(:,1)*T(1,1) + X(:,1) = B(:,1). I.e.,
%             (T(1,1)*S + I)X(:,1) = B(:,1)...

[m,n] = size(B);
I = eye(m,m);
X = zeros(m,n);
X(:,1) = (T(1,1)*S + I)\B(:,1);
for k=2:n
  % Now suppose that we know X(:,1:k-1) for some 1< k < n.
  % Since (by comparing kth columns) 
  %             S*(T(1,k)*X(:,1) + ... + T(k,k)*X(:,k)) + X(:,k) = B(:,k)
  % we see that 
  %                (T(k,k)*S + I)X(:,k) = B(:,k) - S*v 
  %
  % where
  v = X(:,1:k-1)*T(1:k-1,k);
  
  % The product S*v is special since S is upper triangular so we compute
  % B(:,k) - S*v this way...
  for j=1:m
      B(1:j,k) = B(1:j,k) - v(j)*S(1:j,j);
  end
  X(:,k) = (T(k,k)*S + I)\B(:,k);
  % Note that the " + I" is O(m^2), a computation that rivals the 
  % triangular solve flop count. So from the flop point of view should do
  % this via
  %    A = T(k,k)*S; for i=1:m, A(i,i) = A(i,i)+1, end; X(:,k) = A\B(:.k).
  % This was disregarded in grading.
end


% To guaranteed that the solve process doed not break down, it must be the
% case that no diagonal element of S is the negative reciprocal of a
% diagonal entry from T. (This ensures that no diagonal of T(k,k)*S + I
% is zero.
