Taylor expanding $f(y+\epsilon U1 + {\epsilon}^{2} U2,t,\epsilon)$ in $\epsilon$ How would one Taylor expand $\epsilon f(y+\epsilon U1 + {\epsilon}^{2} U2,t,\epsilon)$ in $\epsilon$? 
Somehow the professor obtained the first few terms to be: $\epsilon f(y+\epsilon U1 + {\epsilon}^{2} U2,t,\epsilon)$ $\approx$ $\epsilon$ $f(y,t,0)$ + ${\epsilon}^{2} f_{y}(y)U1 + {\epsilon}^{2} f_{\epsilon}(y,t,0)$? 
 A: The Taylor expansion of a function $g(\epsilon)$ around $\epsilon = 0$ to second order is given by
\begin{equation}
 g(\epsilon) \approx g(0) + \epsilon \, \frac{\text{d} g}{\text{d} \epsilon}(0) + \frac{\epsilon^2}{2}\,\frac{\text{d}^2 g}{\text{d} \epsilon^2}(0). \tag{1}
\end{equation}
The problem here is that the function $g$ depends in a rather complicated way on $\epsilon$, namely
\begin{equation}
 g(\epsilon) = \epsilon\,f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon). \tag{2}
\end{equation}
The first term in $(1)$ is rather easy to calculate: we see that $g(0) = 0 \cdot f(y,t,0) = 0$. Now, for the second term in $(1)$, we need to calculate the derivative of $g$ to $\epsilon$. Using $(2)$, we see that
\begin{equation}
 \frac{\text{d} g}{\text{d} \epsilon} = f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon) + \epsilon \frac{\text{d}}{\text{d}\epsilon} \Big[f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon)\Big] \tag{3}
\end{equation}
by the product rule. The derivative of $f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon)$ to $\epsilon$ is still quite complicated. For the moment, let's denote it as
\begin{equation}
 \frac{\text{d}}{\text{d}\epsilon} \Big[f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon)\Big] =: h(\epsilon).\tag{4}
\end{equation}
Luckily, we don't need to calculate it at this stage: whatever it is, we see that
\begin{equation}
 \frac{\text{d} g}{\text{d} \epsilon}(0) = f(y,t,0) + 0\cdot h(0) = f(y,t,0).
\end{equation}
For the third term in $(1)$ though we do have to get our hands a bit dirty. We need to calculate the second derivative of $g$ to $\epsilon$. This can be done by taking the $\epsilon$-derivative of $(3)$:
\begin{align}
 \frac{\text{d}^2 g}{\text{d} \epsilon^2} = \frac{\text{d}}{\text{d} \epsilon} \Big[\frac{\text{d} g}{\text{d} \epsilon}\Big] &= \frac{\text{d}}{\text{d} \epsilon} \Big[f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon) + \epsilon\, h(\epsilon)\Big] \\ &= \frac{\text{d}}{\text{d} \epsilon} \Big[f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon)\Big] + h(\epsilon) + \epsilon \frac{\text{d} h}{\text{d}\epsilon}\\ &= 2\, h(\epsilon) + \epsilon \frac{\text{d} h}{\text{d}\epsilon}. \tag{5}
\end{align}
When evaluated at $\epsilon = 0$, we obtain
\begin{equation}
\frac{\text{d}^2 g}{\text{d} \epsilon^2}(0) = 2\,h(0) + 0\cdot \frac{\text{d} h}{\text{d} \epsilon}(0) = 2\,h(0).
\end{equation}
Now, what is $h(\epsilon)$? For this, we need the multidimensional chain rule. Namely, we have function of three variables: $f(x_1,x_2,x_3)$. In addition, in the place of these variables, we substitute functions of $\epsilon$. Denoting
\begin{align}
 z_1(\epsilon) &= y + \epsilon U_1 + \epsilon^2 U_2,\\
 z_2(\epsilon) &= t,\\
 z_3(\epsilon) &= \epsilon,
\end{align}
we see that
\begin{equation}
 f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon) = f\big(z_1(\epsilon),z_2(\epsilon),z_3(\epsilon)\big).
\end{equation}
Now we can use the multidimensional chain rule to calculate $h$: we see that
\begin{equation}
h(\epsilon) = \frac{\text{d}}{\text{d}\epsilon} \Big[ f\big(z_1(\epsilon),z_2(\epsilon),z_3(\epsilon)\big)\Big] = \frac{\text{d} z_1}{\text{d}\epsilon}\cdot f_{x_1}(z_1(\epsilon),z_2(\epsilon),z_3(\epsilon)) + \frac{\text{d} z_2}{\text{d}\epsilon}\cdot f_{x_2}(z_1(\epsilon),z_2(\epsilon),z_3(\epsilon)) + \frac{\text{d} z_3}{\text{d}\epsilon}\cdot f_{x_3}(z_1(\epsilon),z_2(\epsilon),z_3(\epsilon)),
\end{equation}
where $f_{x_1} = \frac{\partial f}{\partial x_1}$ is the partial derivative of $f$ to its first variable, $x_1$.
In our case, we have
\begin{align}
 \frac{\text{d}z_1}{\text{d}\epsilon} &= U_1 + 2\epsilon U_2,\\
 \frac{\text{d}z_2}{\text{d}\epsilon} &= 0,\\
 \frac{\text{d}z_3}{\text{d}\epsilon} &= 1,
\end{align}
so
\begin{equation}
 h(\epsilon) = \left(U_1 + 2 \epsilon U_2\right) f_{x_1}(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon) + f_{x_3}(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon),
\end{equation}
hence
\begin{equation}
 h(0) = U_1 f_{x_1}(y,t,0) + f_{x_3}(y,t,0).
\end{equation}
Combining the above expressions, we see that
\begin{equation}
\epsilon\,f(y+\epsilon U_1 + \epsilon^2 U_2,t,\epsilon) \approx 0 + \epsilon\,f(y,t,0) + \epsilon^2 \Big(U_1 f_{x_1}(y,t,0) + f_{x_3}(y,t,0)\Big).
\end{equation}
