function Y = sylv_triangular(TA, LB, D)
% solve a triangular Sylvester equation
% TA*Y - Y*LB = D

[m, n] = size(D);
Y = zeros(m, n);

for j = n:-1:1
    for i = m:-1:1
        num = D(i,j) ...
            - TA(i,i+1:end)*Y(i+1:end,j) ...
            + Y(i,j+1:end)*LB(j+1:end,j);
        den = TA(i,i) - LB(j,j);
        Y(i,j) = num / den;
    end
end