How to expand in taylor series a composition of functions?

I would like to expand in taylor series a function of type : f(x+f(x)) around x=a in the case where f(a)=0.

(%i1) atvalue(f(x),[x=a],0)$

The direct calculus yields :

(%i2) taylor(f(x+f(x)),x,a,2);

(%o2)/T/ f(a)+(at('diff(f(f(x)+x),x,1),x=a))*(x-a)+((at('diff(f(f(x)+x),x,2),x=a))*(x-a)^2)/2+...

If I define a intermediate function :

(%i3)define(tf(x),taylor(f(x),x,a,2))$

Then a expand in Taylor series I get :

(%i4) taylor(f(x+tf(x)),x,a,2);

(%o4) 0+...

I expect the following result : f(1+f'(a))f'(a)(x-a)+(x-a)^2 f''(a)[f'(a)+(1+f'(a))^2/2]+o(x-a)^2

How could I solve this problem?

Solution

You can use gradef to simplify notation.

gradef(f(x),  f1(x)) $
gradef(f1(x), f2(x)) $

atvalue(f(x), x = a,  0) $

e: f(x+f(x)) $
e: taylor(e, x, a, 2) $
e: expand(e, 0, 0)$ /* 'taylor' form to ordinar expression */
e: ev(e, nouns);    /* f(a) to 0 */

returns

          2                                          2
       (f1 (a) f2(a) + 3 f1(a) f2(a) + f2(a)) (x - a)
(%o7) -----------------------------------------------
                              2
                                                          2
                                                     + (f1 (a) + f1(a)) (x - a)