  function X = TriMatSol(S,T,B)
% S is nxn and unit lower triangular
% T is nxn and unit upper triangular
% B is nxn 
% X is nxn and xsatisfies SXT + X = B
%
[n,n] = size(B);
X = zeros(n,n);
% Consider the first column in the equation  SXT + X = B:
%          SX(:,1) + X(:,1) = B(:,1)
% i.e., (S+I)X(:,1) = B(:,1)
Stilde = S + eye(n,n);
X(:,1) = Stilde\B(:,1);
% Now assume k>1 and that X(:,1:k-1) is knowm.
% Let Y = XT and note that the kth column of Y  is given by 
%   Y(:,k) = XT(:,k)
%          = X(:,1:k-1)T(1:k-1,k) + X(:,k)T(k,k)
%          = X(:,1:k-1)T(1:k-1,k) + X(:,k)
% Compare the k-th column of B = SXT + X = SY + X
%   B(:,k) = SY(:,k) + X(:,k) 
%          = SX(:,1:k-1)T(1:k-1,k) + SX(:,k) + X(:,k)
% Solve this for X(:,k)
%   (S + I)X(:,k) = B(:,k) - SX(:,1:k-1)T(1:k-1,k)
for k=2:n
    X(:,k) = Stilde\(B(:,k) - S*(X(:,1:k-1)*T(1:k-1,k)));
end

% -1 if X(:,1:k-1)*T(1:k-1,k) is computed as a sum of vectors.

