cf_matrixfunction(fct,A,b,n,varargin)

function [v] = cf_matrixfunction(fct,A,b,n,varargin)
% compute the product v = f(A)*b without the explicit computation of f(A)
% The matrix A has to be negative semidefinite (that is, all eigenvalues
% should be negative)

% First the method computes a Caratheodory-Fejer approximation, that is, a
% a uniform rational approximation of the function f of type (n,n) on 
% the negative real axis. 

% The matrix A has to be real and n should be even
% The vector b has to be real, too.

[zi,ci,r_inf] = cf_realaxis(fct,n,varargin{:});

% identity matrix
I = speye(size(A));
v = real(r_inf*b);
for i=1:n/2
    v = v + 2*real(ci(i)*((A-zi(i)*I)\b));
end