funm_randomized

function F = funm_randomized(A,fun)
%FUNM_RANDOMIZED  Evaluate general matrix function using
% randomized approximate diagonalization method of Davies (2007).
tol = 8*eps(A(1,1));  % Tolerance for single or double precision A.
E = randn(size(A)); 
[V,D] = eig(A + (tol*norm(A,'fro')/norm(E,'fro'))*E);
F = V*diag(fun(diag(D)))/V;

The code is easy to convert to work with any precision (not only single or double).
The only one change is to create random matrix of the same type as A:

E = randn(size(A),class(A));

Now our code has became precision-independent and can be used with double precision:

>> A = rand(5,'double');
>> norm(funm_randomized(A,@exp)-expm(A),1)/norm(A,1)
ans =
      7.75329562589415e-15

or quadruple precision:

>> mp.Digits(34);
>> A = rand(5,'mp');
>> norm(funm_randomized(A,@exp)-expm(A),1)/norm(A,1)
ans = 
    7.578684267133064009245947356284363e-33

or any other:

>> mp.Digits(1000);
>> A = rand(5,'mp');
>> norm(funm_randomized(A,@exp)-expm(A),1)/norm(A,1)
ans = 
    3.16466380.....10447022701e-998

Thanks, Pavel - that's an excellent point. E should indeed be of the same precision as A and your change achieves that very neatly.