phipade - Evaluate phi functions using (6,6)-Pad� approximations.
 
  SYNOPSIS:
     phi_k                     = phipade(z, k);
    [phi_1, phi_2, ..., phi_k] = phipade(z, k);
 
  DESCRIPTION:
    This function evaluates phi functions needed in exponential
    integrators using diagonal Pad� approximants (currently (6,6)-Pad�).
    We define the phi functions according to the integral representation
 
       \phi_k(z) = \frac{1}{(k - 1)!} \int_0^1 e^{z (1-x)} x^{k-1} dx
 
    for k=1, 2, ...
 
  PARAMETERS:
    z - Evaluation point.  Assumed to be one of
          - 1D vector, treated as the main diagonal of a diagonal matrix
          - sparse diagonal matrix
          - full or sparse matrix
    k - Which phi function(s) to evaluate.
        Index (integer) of the (highest) phi function needed.
 
  RETURNS:
     phi_k                     =      \phi_k(z)
    [phi_1, phi_2, ..., phi_k] = DEAL(\phi_1(z), \phi_2(z), ..., \phi_k(z))
 
  NOTES:
    When computing more than one phi function, it is the caller's
    responsibility to provide enough output arguments to hold all of the
    \phi_k function values.
 
    For efficiency reasons, phipade caches recently computed function
    values.  The caching behaviour is contingent on the WANTCACHE
    function and may be toggled on or off as needed.