clear i
N = 10;
A = rand(N)+i*rand(N);
expmA = expm(A);
I = eye(N);
disp('phi1m')
ref = A\(expmA-I);
norm(phi1m(A)-ref)/norm(ref)
disp('phi2m')
A2 = A^2;
ref = A2\(expmA-I-A);
norm(phi2m(A)-ref)/norm(ref)
disp('phi3m')
A3 = A2*A;
ref = A3\(expmA-I-A-A2/2);
norm(phi3m(A)-ref)/norm(ref)
disp('phi4m')
ref = (A2*A^2)\(expmA-I-A-A2/2-A3/6);
norm(phi4m(A)-ref)/norm(ref)