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

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