function x = ProdTriSol(S,T,lambda,b)
% S and T are n-by- and upper triangular.
% b is n-by-1.
% lambda is a scalar and ST - lambda*I is nonsingular.
% x solves the nonsingular system (ST - lambda*I)x = b.

n = length(b);
x = zeros(n,1);
w = zeros(n,1);
x(n) = b(n)/(S(n,n)*T(n,n)-lambda);
w(n) = T(n,n)*x(n);
for k= n-1:-1:1
% Generate x(k:n) and w(k:n) so
% (S(k:n,k:n)T(k:n,k:n) - lambda I)x(k:n) = b(k:n)
% w(k:n) = T(k:n,k:n)x(k:n)
mu = T(k,k+1:n)*x(k+1:n);
x(k) = (b(k) - S(k,k)*mu - S(k,k+1:n)*w(k+1:n))/...
(S(k,k)*T(k,k)-lambda);
w(k) = T(k,k)*x(k) + mu;
% A pair of length n-k dot products per pass thru the loop
% Notice that if you are careless and use

%x(k) = (b(k) - S(k,k)*T(k,k+1:n)*x(k+1:n) - S(k,k+1:n)*w(k+1:n))/...
% (S(k,k)*T(k,k)-lambda);
%w(k) = T(k,k)*x(k) + T(k,k+1:n)*x(k+1:n);

% then you have an extra length n-k vector scale (-1) and an extra length
% n-k dot product (-1). (Note: scalar*vector'*vector is an order of
% magnitude more work than scalar*(vector'*vector)

end