function X = funm_parlett(f, A)
[Q, U] = schur(A, 'complex');

n = size(U, 1);
S = zeros(n, n);
for j = 1:n
    S(j,j) = f(U(j,j));
    for i = j-1:-1:1
        S(i,j) = ((U(i,i+1:j)*S(i+1:j,j) - S(i,i:j-1)*U(i:j-1,j))) / ...
            (U(j,j) - U(i,i));
    end
end
X = Q * S * Q';