Diff and Taylor to return expressions in terms of partial derivatives?

For the derivative, maybe you would like to use pdiff (positional derivative package)

  (%i1) load(pdiff)$

  (%i4) subst(tau = epsilon*t, f0(t, tau) + epsilon*f1(t, tau) + epsilon^2*f2(t,tau));
  (%o4) epsilon^2*f2(t,epsilon*t)+epsilon*f1(t,epsilon*t)+f0(t,epsilon*t)

  (%i5) diff(%,t);
  (%o5) epsilon^2*(f2[(1,0)](t,epsilon*t)+epsilon*f2[(0,1)](t,epsilon*t))+epsilon*(f1[(1,0)](t,epsilon*t)+epsilon*f1[(0,1)](t,epsilon*t))+f0[(1,0)](t,epsilon*t)+epsilon*f0[(0,1)](t,epsilon*t)

I think pdiff doesn't allow subscripted function names.

--Barton
________________________________


(1) Introduce a slow variable: tau=epsilon*t

(2) Introduce the assumption that: f(t, tau, epsilon) = f[0](t, tau) + epsilon*f[1](t, tau) + epsilon^2*f[2](t,tau) + ...
In the textbook f[i] --> f superscript i,
This is like a Taylor series expansion, except that the variable t is never multiplied on the outside to avoid secular terms

(3) df[0]/dt = partial(f[0], t) * partial(t,t) + partial(f[0], tau) * partial(tau, t) = f[0]_1 + epsilon*f[0]_2
Where the underscore denotes partial with respect to the ith argument

(4) Thus
df/dt = f[0]_1 + epsilon*(f[0]_2 + f[1]_1) + epsilon^2*(f[1]_2 + f[2]_1) + epsilon^3*(f[2]_2 + f[3]_1) + ...
d^2f/dt^2 = f[0]_11 + epsilon*(f[0]_12 + f[0]_21 + f[1]_11) + epsilon^2*(f[0]_22 + f[1]_12 + f[1]_21 + f[2]_11) + ...

So the question. How can I tell Maxima to replicate the above results easily in Maxima?