phipade - Evaluate phi functions using (6,6)-Pad� approximations.

SYNOPSIS:
[phi_1, phi_2, ..., phi_k] = phipade(z, k);

DESCRIPTION:
This function evaluates phi functions needed in exponential
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.