function Y = sylv_triangular(UA, LB, E)

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

for j = n:-1:1
    for i = m:-1:1
        Y(i,j) = E(i,j) / (UA(i,i) - LB(j,j));
        E(1:i-1,j) = E(1:i-1,j) - UA(1:i-1,i) * Y(i,j);
        E(i,1:j-1) = E(i,1:j-1) + Y(i,j) * LB(j, 1:j-1);
    end
end