Composition of Taylor Series Suppose I have smooth functions $f,g,y_0$ and $y_1$ from $\mathbb{R}$ to $\mathbb{R}$, such that
$$y_1(x) = y_0(x) - \epsilon g(y_0(x))$$
Then I consider
$$f(y_0(x)) = f(y_1(x) + \epsilon g(y_0(x)))$$
Is there a closed form expression for the Taylor series in the small parameter $\epsilon$ in terms of derivatives of $f$ and $g$ and only the function $y_1$? 
The first few terms are
$$f(y_0) = f(y_1) + \epsilon f'(y_1) g(y_0) +  \frac{1}{2}\epsilon^2 f''(y_1)g^2(y_0) +..$$
Where we interpret $f(y_0)$ as $f|_{y_0(x)}$ and treat $x$ fixed. Then we can again replace the $y_0$ in $g(y_0)$ with
$$g(y_0) = g(y_1)+ \epsilon g'(y_1)g(y_0) +...$$
giving
$$= f(y_1) + \epsilon f'(y_1) [g(y_1) + \epsilon g'(y_1)g(y_0) + ... ]$$
$$+  \frac{1}{2}\epsilon^2 f''(y_1)[g(y_1) + \epsilon g'(y_1)g(y_0) + ... ]^2 +...$$
Continuing to replace the $y_0$ with $g(y_0)$ like this and grouping terms gives
$$f(y_0) = f(y_1)+ \epsilon [f'g](y_1) + \epsilon^2[f'g'g + \frac{1}{2}f''g^2](y_1) + \epsilon^3[f'g'^2g + \frac{1}{2}f'g''g + \frac{1}{2}f''g'g + \frac{1}{6}f'''g^3](y_1) + O(\epsilon^4)$$
But is there some way to write this as a more compact sum like
$$f(y_0) \sim f(y_1) + g(y_1)\sum_{n=1}^\infty\sum_{m=0}^n \epsilon^n \alpha(n,m)f^{(n)}g^{(n-m)}(y_1)$$
I am having trouble identifying the pattern. I know there will be some product involved as well.
Edit:
Thinking about it some more it may suffice to just set $f=id$ and consider
$$y_0 = y_1 + \epsilon g(y_0)$$
$$y_0 = y_1 + \epsilon g(y_1 + \epsilon g(y_0))$$
$$y_0 = y_1 + \epsilon g(y_1 + \epsilon g(y_1 + \epsilon g(y_0)))$$
$$y_0 = y_1 + \epsilon g(y_1 + \epsilon g(y_1 + \epsilon g(y_1 + ...)))$$
and somehow use the chain rule
$$[f_1\circ f_2 \circ .... \circ f_n]' = \prod_{i=1}^n(f'_{i}\circ f_{i+1}\circ ...\circ f_n)$$

It looks like there is a common factor of εg(d/dy) being used as both a function and a multiplicand. 
It's being used as a function through composition -repeatedly taking the derivative and multiplying by g N times-, and as a multiplicand -taking the derivative N times then multiplying by g to the N.
Each εⁿ term looks to be made up of the difference between these two operations. Not sure of the factor that each is scaled by before taking the difference. (Repeatedly taking the derivative of each term and then multiplying by εg to obtain the next term is off by a constant factor of N and an additional amount of (d^n/dy^n)f*(g^n) scaled by N!.)
