function E = phi1m(A)
%PHIM1  Matrix phi function via Pade approximation.
%   M. Hockbruck, C. Lubich and H. Selhofer,
%   Exponential integrators for large systems of differential equations,
%   SIAM J. Sci. Comp. (1996)
%
% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% Pade approximation for phi(A)
ID = eye(size(A));
X = A;
q=6;
c=[1/26 5/156 1/858 1/5720 1/205920 1/8648640];
E = ID + c(1)*A;
D = ID - 2*q*c(1)*A;
p = 1;
for k = 2:q
%   c = c * (q-k+1) / (k*(2*q-k+1))
   X = A*X;
   cX = c(k)*X;
   E = E + cX;
   if p
     D = D + (k+1)*cX;
   else
     D = D - (2*q-k+1)*cX;
   end
   p = ~p;
end
E = D\E;

% Undo scaling by repeated squaring
EXP = A * E + ID;
for k=1:s
  E = (EXP + ID) * E / 2;
  EXP = EXP * EXP;
end
